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ABSTRACT 


The behavior of the electric field together with the electron and ion densities in the 
vicinity of a non-emitting, plane anode 1s investigated. The selected approach involves non- 
linear analysis techniques on the continuum equations for steady-state, isothermal conditions 
where both ionization and two-body recombination are included. Ions, created through 
electron bombardment of neutral atoms, are repelled toward two stagnation regions: within 
or near the sheath boundary and near the plasma interface. These equilibria form as a result 
of the chemistry present: recombination establishes the latter while ionization stipulates the 
former. As presented, the sheath is fundamentally unstable - ions are driven toward the 
negative electrode. Using nitrogen data for a numeric example, the following observations 
are made: a sufficiently strong applied electric field pushes the ion density toward that of the 
electrons through a well - a constrictive phenomenon. Both a transition region, dominated 
by density gradients, and a diffusion-driven zone are found to move the system toward the 
plasma interface. The characteristics of this process are influenced by the applied electric 
field, but the instability of the chemistry-induced stagnation regions precludes numeric 
convergence. Insufficient dissipation may prevent the stability of the anode fall model as 
presented. Suggested improvements to the model descriptions include considering the effects 
of temperature gradients, magnetic fields, three-body recombination, diffusion written in 


terms of the electric field, multi-dimensionality and/or time-dependencies. 
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I. INTRODUCTION 


Whenever an object starts moving in a given direction, a second object must be 
present which accepts the reaction of the force that accelerates the first, so says Newton's 
third law of motion. When a spacecraft moves through empty space, the accelerating forces 
act between the vehicle and the mass of the propellant expelled. In 1895, the distinguished 
Russian rocket pioneer Konstantin Eduarovitch Tsiolkovskii derived the famous *Tsiolkovskii 
equation," the basis of all theoretical work on rocket propulsion. Relating the burnout 
velocity of a rocket with its exhaust velocity and its mass ratio, this “rocket equation” shows 
how important it is for the exhaust velocity of a rocket motor to be as high as possible. [Ref. 
1: pp. 2-3] 

Since mankind”s advent into space, several types of space flight propulsion systems 
have been developed to include chemical, nuclear, electric and solar propulsion devices with 
the majority of spacecraft to date employing chemical rockets. While the science and 
technology of chemical rockets developed at a rapid pace, the idea of using the benefits of 
Maxwell's relations to eject charged particles from a thrust chamber progressed slowly with 
long intervals between periods of creative thought. In 1945, Herbert Radd recommended the 
use of fast electrically accelerated particles to reduce fuel mass and in 1946, Jakob Ackeret 
analyzed the performance capability of a spacecraft by arguing that the total kinetic energy 
carried away by the exhaust beam corresponded to a reduction of the mass of the vehicle 
according to E=mc*. Hence, assuming that a certain fraction of the vehicle mass is converted 
into kinetic beam energy, the total propellant mass can be calculated which in turn leads to 
a maximum terminal velocity of the spacecraft, a result which is independent of particle size, 
time of propulsion, conversion method and magnitude of the vehicle. In essence, the 
optimum exhaust velocity is determined simply by maximizing the terminal velocity. [Ref. 1: 
pp. 2-3] 

Propulsive devices that use Maxwell’s Law’s differ from chemical systems: the latter 
is characterized by having the same energy and propulsive sources, while in the former the 


energy conversion method differs from that of propulsion. In essence, electric propulsion 


systems allow greater carrying capacity with shorter flight times. For example, the Hohmann 
transfer time from a Low-Earth-Orbit (LEO) to the orbit of Pluto requires approximately 45 
years if chemical means are considered; an electrically propelled vehicle, in comparison, 
requires approximately three years. [Ref. 1: p. viii]. The large difference in mass payload ratio 
(final mass/initial launch mass) obtained 1s highlighted by an analysis for a Mars mission: a 
chemical system using a high impulse Hohmann trajectory from LEO delivers approximately 
1095 to 1896 of the launch mass to the Martian surface [Ref. 2: pp. 152], while an electric 
system using a low impulse spiral trajectory could deliver 20% to 60% of the launch mass, 
depending on the desired transit time. With effective exhaust velocities as high as 10,000 m/s, 
electric propulsion thus offers the performance envelope needed for manned interplanetary 
missions. [Ref . 2: pp. 37] 

Electric propulsion systems are divided into three categories: 1) electrostatic, where 
acceleration is achieved through the interaction of coulombic fields with charged propellant 
particles, 2) electrothermal, where propellant heated electrically is expanded 
thermodynamically through a nozzle, and 3) electromagnetic, where acceleration is achieved 
through the interaction of the Lorentz force with highly ionized gases (plasma). Of the types 
of electric propulsive thrusters available, electrothermal devices include the resistojet and 
arcjet. It should be noted that the current necessary to form the plasma in these devices also 
allows the magnetic field to accelerate the gas. These systems (and for that matter, any 
device using a thermal nozzle) are, however, limited by the amount of heat that can be added 
to the flow and are thus ultimately dependent on the materials used. As electromagnetic 
(EM) thrusters rely on Lorentz interactions to accelerate the flow, these devices are not 
limited in the sense of electrothermal systems and therefore offer different operational 
characteristics. The various types of EM thrusters available include the ion engine and plasma 
thrusters (wide variety to include electrodynamic, magnetoplasma-dynamic or MPD, pulsed- 
plasma, Hall accelerator, Lorentz force accelerator, pulsed coaxial thruster, plasma arcjet, 
pulsed ablative thrusters, to name just a few). [Ref. 2: pp. 565-596] The MPD thruster is 
relevant to this work. 


MPD systems are a relatively recent addition to the electric propulsion family, 


originally discovered when peculiar behavior of an electrothermal arcjet thruster was 
investigated. Researchers found that without significant mass flow-rate applied, the arcjet 
continued to operate in a mode clearly unique: pressure within the thrust chamber was too 
low for the device to be choked, and so it was theorized that an induced magnetic field from 
the high arc current accelerated the plasma flow - the magnetoplasmadynamic thruster was 
“born.” [Ref. 3] 

The area of research to which MPD thrusters are applied is termed 
magnetohydrodynamics (MHD), the field of physics that investigates how moving conductive 
(electrolytic) fluids behave while subjected to applied electric or magnetic fields. The fluids 
involved must be conductive, or ionized sufficiently to pass a current, and in the case of an 
MPD, the fluid must reach a plasma state. MHD devices behave analogously to motors or 
generators and can either add or remove energy from the fluid in exchange for electrical 
energy; systems developed include the MHD tunnel drive for submarines and the MHD 
turbine. Magnetoplasmadynamic (MPD) thrusters are currently being considered for 
application to primary propulsion system on robotic and piloted interplanetary missions; the 
high specific impulses offered by MPD devices can reduce launch mass requirements by over 
a factor of two compared to chemical system if the power to thrust conversion efficiency is 
over 50%. [Ref. 4] 

MPD thrusters are classified by considering two major features: 1) the power supply 
which defines the continuous, quasi-steady and pulsed MPD systems, and 2) the magnetic 
field which characterizes the self-field and applied-field MPD thrusters [Ref. 3]. A 
continuous, or steady-state, thruster is one in which the propellant flow and current 
throughout the device are continuous, thereby allowing the system to reach a high 
temperature equilibrium state. The quasi-steady thruster does not reach equilibria since it 
uses a continuous gas flow with rapidly pulsed current, through which high amounts of energy 
are provided in short bursts. The pulsed thruster is similar to the quasi-steady MPD system 
except it uses less frequent pulses; equilibrium flow is also not reached and consequently low 
temperature materials may be used for constructing such a device. The applied- and self-field 


thrusters define the way in which the primary accelerating field (magnetic) is affected. The 


applied-field system employs a large current in an external coil, or large permanent magnets, 
to accelerate the plasma whereas in the self-field device the current not only ionizes the gas 
but also creates the field that accelerates the resulting plasma. 

The magnetoplasmadynamic thruster was originally discovered in the 1960's, and yet 
much work remains to understand its fundamental operation. For example, the steady state, 
self-field MPD thruster operated by the University of Southampton, U.K., was originally 
developed to simulate the atmosphere effects on structures in LEO, an environment 
dominated by highly reactive atomic oxygen. The MPD is suited to simulate such conditions 
as it offers a unique combination of high thrust and high exit velocity. The Southampton 
MPD presently operates at power levels in the range 6-17 kW but suffers from electrode 
erosion, instabilities, a contaminated beam and short run-times. Work on this thruster has 
identified electrode erosion as a major problem, making the device unsuitable for its intended 
function. Specifically, the mode of erosion differs from that seen at other institutions and the 
levels of erosion are higher and exhibit a rate increase not seen with other devices. [Ref. 5] 
Electrode erosion, current spotting, frozen flow loses and electrode power deposition are just 
several factors that limit the performance of MHD devices [Ref. 6, 7]. Specifically, anode 
power deposition is the single largest power loss mechanism in MPD thrusters operating at 
sub-megawatt power levels. [Ref. 8, 9, 10, 11, 12] 

To gain further insight to MPD phenomena, computer codes that accurately describe 
observed data from steady-state MPD thrusters have been developed [Ref. 13, 14, 15]. 
However, these codes do not adequately describe observed data from quasi-steady thruster 
experiments. It has been suggested that the lack of proper electrode modeling (i.e., sheaths 
and fall potentials) in these codes may explain this discrepancy [Ref. 16]; perhaps the difficult 
set of coupled, nonlinear partial differential equations involved has limited analytical work to 
model the anode sheath and ambipolar regions. For example, Hugel [Ref. 17] and 
Subramanian [Ref. 18] address the influence of the sheath region, but do not model the 
electric field, temperature or sheath fall voltage. Hence, the focus of this work: analysis of 


the anode voltage drop region. 


Il. LITERATURE REVIEW 


The principal voltage loss mechanism near an electrode can be divided into two main 
classes: ohmic and sheath losses. Ohmic resistivity occurs because of the finite conductivity 
of a real plasma, particularly thermal boundary layers, degree and kinetics of ionization and 
Joule heating. Sheath drops are caused by Debye shielding which forms a space charge field 
adjacent to the electrode [Ref. 19: pp. 200-204]. Furthermore, material problems restrict the 
temperatures at which the electrodes can operate. In many cases cooling of the electrodes 
Is required since the plasma, in order to maintain a high ionization, must be hotter than the 
working temperatures of most materials. This temperature difference between the electrodes 
and the plasma further aggravates the voltage losses because of the presence of the thermal 
boundary layers. As has been shown, voltage losses can exceed 5096 of the total power 
output, with sheath drops accounting for a large fraction of this loss. [Ref. 6, 7] 

The tool most commonly used for measurement of local plasma properties is the 
electrostatic probe, a small electrode which may be biased positively or negatively with 
respect to the plasma in which it is immersed. Such work is relevant to non-emitting anodes 
which, though heavily biased, also probe and disturb the plasma locally. As a result, most 
analytical work on sheath phenomena is embodied in probe theory investigations where the 
effects of an electrode (probe) locally disturb a quiescent plasma. [Ref. 20,21: pp. 1-7, pp. 
8-9] 

In 1924 Langmuir pioneered the use of electric probes, giving the first sound 
mathematical basis to this diagnostic technique. The analysis of electrostatic probes in 
collision-dominated (higher pressure) plasmas began with Davydov and Zmanovskaja [Ref. 
22] who assumed that a quasi-neutral diffusion controlled region extended to within one mean 
free path of the probe where it was matched to the edge of a free-fall sheath; a transition 
region between the two regimes was not attempted, however. Ecker, et al. [Ref. 23] 
attempted matching the solution in the collision-dominated region directly to that in the free- 
fall sheath using variational principles, work which was criticized by Cohen [Ref. 24] based 


on mathematical grounds: something was forced to fit via some artificially imposed boundary; 


matching between separate regions was accomplished through vanational principles and 
consequently was not asymptotic in nature [Ref. 24]. 

Existing sheath solutions are limited to special case simplifications. Several solutions 
are available for spherical probes by Kiel [Ref. 25], whose analysis neglected diffusion, Barad 
and Cohen [Ref. 26], Su and Lam [Ref. 27], as well as by Cohen [Ref. 24, 28]. The latter 
group carried out the first completely systematic analyses of spherical probes in collision- 
dominated plasmas. These researchers assumed the continuum-equations valid everywhere 
in the plasma right up to the probe surface; their results showed both the sheath and the quasi- 
neutral regions appeared naturally from the diffusion equations. In these efforts, two major 
assumptions prevail: the plasma was assumed so lightly ionized that only charge-neutral 
interactions were accounted, and the continuum equations were postulated with an isothermal 
formulation. For cool, nonequilibrium plasmas, however, charge-charge collisions must be 
considered [Ref. 29: p. 165]. To that end, McKee and Mitchner [Ref. 30] investigated a 
collisionless sheath with constant transport properties that also included ionization as well as 
recombination effects in the ambipolar region. Exploring spherical probe theory further, 
Cohen and Schweitzer [Ref. 31] then extended Cohen's original asymptotic analysis [Ref. 24] 
for the entire region, an analysis which allowed for weak ionization and recombination effects 
while also assuming constant transport coefficients. Su and Sonin [Ref. 32] next formulated 
a spherical probe theory by considering the effects of electron-ion collisions by approximating 
diffusion coefficients as a function of the neutral particle density, a density assumed as 
constant. Barad and Cohen [Ref. 26] then presented a continuum theory (for spherical 
probes) in which non-constant transport coefficients were used, coefficients based on charge- 
charge collisional assumptions. This analysis, however, did not incorporate ionization or 
recombination effects. In work similar to Lam's original theory for flow of weakly ionized 
gases [Ref. 28], Stahl and Su [Ref. 33] use the same approach of separating the sheath, 
ambipolar region and free stream; their methods prove the existence of a sheath on a flat 
probe, yet plasma chemistry is not considered (ionization and recombination). It should be 
noted that both Su and Lam [Ref. 27], for spherical probes, and Stahl and Su [Ref. 33], for 


flat probes, consider transition regions to match the inner sheath to the outer quasi-neutral 


or ambipolar zones - albeit both research groups neglected plasma reactivity in their 
formulations. In related work, Godyak and Sternberg [Ref. 34] applied the ionization 
potential as the boundary condition for both the plasma and the sheath in an effort to 
understand the electrodynamic properties of the plasma-sheath boundary for the cathode; their 
analysis, however, assumed constant ionization for the plasma reactivity with recombination 
neglected. 

It may be inferred from these efforts that anode-voltage-drop-region analyses remain 
insoluble due not only to the nonlinear nature of the equations, but also to the transitional 
nature of particle motion in most sheaths. In an extensive review of probe theory, Chung, 
Talbot and Touryan [Ref. 20] state that a general solution for determining charge density and 
species temperature for probes small relative to boundary layer thicknesses is not available. 
Additionally, Nasser [Ref. 35] discusses an elementary theoretical approach to the glow 
discharge problem by introducing a set of one-dimensional governing equations (continuum) - 
the electron and ion current equations, the net production of electrons and ions, and the 
Poisson equation. According to Nasser, these equations should have a solution, but most 
attempts to achieve such have failed, with the difficulty lying in the boundary conditions [Ref. 
35: p. 412]. Voltage losses at electrode boundanes and surface erosion, together with sheath 
effects which influence ionization, play an important role in plasma devices and must be 
controlled in designs of practical interest. Particularly, thermal arcjets and MPD accelerators 
deposit between 15% and 80% of the input power into the anode. This presents not only a 
severe performance penalty, but also introduces thermal design problems since the heat thus 
generated must be radiated away from the thruster surfaces [Ref. 7, 9]. 

Completely satisfactory solutions to the anode voltage drop phenomena do not exist 
for several plasma conditions; one such case involves steady, collisional, low-temperature 
isothermal discharges [Ref. 36]. For instance, the effect of several parameters on the region 
require further study, including the influence of temperature, the role of the problem 
dimensionality as well as the effect of plasma chemistry. It remains unclear, moreover, as to 
the exact cause of the anode spots which have been observed for various MPD thruster 


operating conditions, a behavior which appears roughly as periodic oscillations [Ref. 6, 7, 9]. 


A one-dimensional, collisional, equilibrium solution can satisfactorily reproduce the observed 
electric field and charge density distributions for the entire anode voltage drop region where 
the electron temperature equals that of the heavy species [Ref. 36]. However, this model 
cannot describe any decrease in current density away from the surface, or current constriction, 
at the anode surface which might be necessary in non-equilibrium. A two-dimensional model, 
developed by Biblarz and Dolson [Ref. 6], represents these phenomena, in the absence of 
species ionization and recombination, and predicts the voltage drop in the region: the sheath 
accounted for a majority of the anode voltage drop [Refs. 6, 7, 36]. 

Since the pioneering work of Tonks and Langmuir (relating to low-pressure 
discharges), the anode description has remained as one of the classic problems in plasma 
physics. In this work, the high density problem is generalized to any electrode surface in 
contact with a partially-ionized plasma when the current flow is sufficiently small: the 
problem of a non-emitting electrode (probe) in contact with a thermally generated plasma. 
An anode is explicitly considered along with a non-emitting cathode and the governing one- 
dimensional continuum equations are formulated. The electrodynamic behavior of plasma 
reactivity in the electrode regions is then discussed with considerations given to three distinct 
zones of the anode voltage drop regime: inner or space-charged sheath, transition and outer 


or quasi-neutral ambipolar zones. 


HL GOVERNING EQUATIONS AND NUMERICAL DATA 


A. GOVERNING EQUATIONS 

A flat electrostatic probe immersed in a high-density, flowing, partially-ionized gas is 
considered. For plasmas of moderately high density and low degree of ionization (as applied 
to discharge lasers and MHD devices), the plasma region affected by an electrode, 1.e., the 
electrode region, can be considered of the boundary-layer type with the transport of charges 
approximated as one-dimensional. Traditionally, the electrode region is further separated into 
two distinct zones, as delineated by the species densities: a sheath and an ambipolar region, 
both which may be of the order of tenths or even hundredths of a millimeter in extent [Ref. 
6] where the sheath extent is driven by variation of the electric field and the quasi-neutral 
region extent is influenced by the influences of chemistry. 

The problem is formulated by considering the hydrodynamic continuum equations for 
two species, ions and electrons. Induced and applied magnetic fields are neglected with 
energy conservation implicitly satisfied since heat transfer, radiation, Joule heating and surface 
emission effects are at first neglected (for example, B=0, VI =0, VT=0, J-E=0). 

The equations below are assumed valid throughout the entire domain of interest, to 


include the electrode (probe) surface. Table (1) lists the nomenclature used: 


o(Mn,v,) | 
ET ncm ee Cope ee qme Mn VIE. Equation (12a) 
Q (mn, v.) 
mure Ds v,Vv,--qn,E-kT Vn,-mn,f,v,-F, Equation (1b) 
q 
V-E= — (n,- 
e fue Equation (2) 
Vigor (il, ) Equation (3a) 
Vj,=aqv(n,v.) Equation (3b) 
NV ACG Vi ee con.n: Equation (4a) 
Vin wn e non. Equation (4b) 
Vy-Qo => EM. Equation (5) 


Equations (1a) and (1b) describe the momentum-transfer between the ions and electrons with 
neutrals: the left-hand side represents the unsteady and convective contributions, 
respectively, while the right-hand side symbolizes coulombic, pressure, collisional and 
frictional forces, respectively. Next, the electrode is considered as a positive, non-emitting 
surface (anode). Through this process, a flow of negative charges is induced toward the plate 
and a flow of positive charges is repelled away from the plate; the anode, in completing the 
arc current with the cathode, collects incoming electrons while driving all positive charges 
toward the negative surface. This separation of charge in the anode fall region produces a 
divergence of the electric field and a variation of potential away from the plate, as given by 
equation (2) (the variation of potential is readily deduced from Gauss’s equation (2)). The 
current densities are denoted by equations (3a) and (3b) in terms of the species flux of 
charged particles per unit area. Equations (4a) and (4b) are then the continuity equations for 
the ion and electron flux as caused by the ionization and electron-ion (two-body) 
recombination processes: positive ions are produced within the anode fall region by electron 
impact with neutrals and conversely, positive ions combine with negative electrons to form 
neutrals. Equation (5) is the conservation of total charge; the sum of the species currents 
remains constant for a one-dimensional formulation. Hence, from set (1-5), the anode model 


presented in this work is developed with the following assumptions: 


1. Steady state; free-stream plasma is neutral with uniform charge densities (Saha 
equilibrium). 


2. Mean free path of the charged particles is small compared to the sheath thickness 
which in turn is much smaller than the dimensions of the body; the sheath is treated 
as a continuum. - 

3. Chemical equilibrium in the boundary layer, the two-body recombination 
coefficient is held constant over the entire anode fall region as a first approximation 
[Ref. 29: p. 238]. Further, the ionization coefficient is treated as Townsend’s first 
coefficient. 


4. There is no applied magnetic field and induced magnetic field effects are neglected 
(low Magnetic Reynolds number) - consequently no ion slip. 


5. No continuum radiation losses, no ion emission, no significant Joule heating. 


6. Isothermal: VT.-VT;-VT,-0 and T,-T7T,. The true electron temperature will 
either be close to the anode surface temperature or the free stream temperature since 
thermal equilibrium is assumed in the boundary layer. Because of the relatively small 
number of collisions that the electrons suffer in the boundary layer, it 1s frequently 
argued that the electron temperature is close to the free stream temperature [Ref. 37]. 


7. Diffusion velocities of the ions and electrons due to electric field interactions are 
small in comparison with their thermal velocities; diffusion coefficients are assumed 
constant. 


8. End effects, in the case of a planar surface, are neglected. 


9. Under the assumption of a highly biased probe, it 1s possible to ignore convection 
not only in the very thin transition region, but also in the sheath region, even for 
sheaths with the same order of thickness as the boundary layer [Ref. 33]. Assuming 
perfectly elastic collision between the species and neutrals, essentially a constant 
collisional frequency is assumed sufficiently large so that the convective derivative in 
(1) can be neglected. 


10. With the assumptions of steady state and negligible convection, a spatially 
homogenous fluid results. 


For a semi-infinite plane with coordinate outward from the planar positive surface so 


that V=d/dy, set (1-5) thus becomes: 
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Here, equations (3a) and (3b) have been incorporated with equations (4a) and (4b). Using 








| > "A | KT. m». 
species mobilities, i,,- —2— and the diffusion coefficients, D,,- —- (the Einstein 


¡el ie i,e i,e 





relation u,,= = ), species momentum conservation can be re-written in terms of species 


i,e 


flux, I’, ., where Fick’s law of diffusion is a special case (E=0 so that p;.=0). Specifically, (1) 


and (2) become 


|] mw -tu. n, ,E-D, ,Vn, , Equation (6) 


Essentially, diffusion is a random process in which a net flux from dense regions to less dense 
regions occurs simply because more particles begin their motion in the more dense area - the 
flux is proportional to the density gradient. [Ref. 38: pp. 135-138] But as the species travel 
toward the electrode, reactivity occurs - ions and electrons recombine while electron collision 
with neutrals breaks neutrality adding charges to the electrode region; flow neutrality 1s 
broken and the fluid 1s no longer a plasma. Thus, under constant-property conditions, the 
species continuity equations together with Gauss’s-equation describe the flux of charges in 
the vicinity of an electrode or plasma probe. [Ref. 33, 39] With the application of the Einstein 
relation, set (1-5) are re-cast to give the model used in this work, the geometry for which is 


shown in Figure (1): 
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Figure 1: Geometry of Anode Fall Region 
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To obtain the mathematical model used in this work, the definition of the diffusion coefficient 
is used to re-write the species elastic collision frequencies, £=[kT,]/[D,m,] for isothermal 
conditions, and the current density is defined as j=qn,v, where s=1, e. The resulting equations 


(1-5) are then solved for their highest derivative terms: 
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In the formulation of equations (7-9), 14 variables result (three universal constants 
included). Before analysis of this system, the equations must be written in terms of 
dimensionless variables and proper parameters. For this work, both the Buckingham Pi 
Theorem and Fractional Analysis techniques were applied (Appendix A). The Buckingham 
method, although a useful tool for revealing the dimensionality of the problem, does not bring 
forth the physical effects which influence system behavior. On the other hand, fractional 
analysis coupled with the use of the governing equations, begins to give physical meaning to 
the dimensionless parameters obtained. 

Using undisturbed plasma values for normalizing the variables gives three primary 


length scales for the system modeled by equations (7-9) since 


eam cam os Jo. So Jo for the undisturbed, neutral 


plasma: L,, L, =L, =L,, L, a =L, . The results are summarızed in Table 2 with 
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the length scales listed in Table 3: 


Governing Equations Non-dimensional Non-dimensional Parameters 
[Equations (7-9)] Parameters via via Undisturbed Plasma 
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Sheath: Transition Zone: 
Space-Charge Dominant Number Density Gradient 
Dominant 





Table 3: Length Scales via Fractional Analysis 
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It should be noted that the length scales Lg, L, and L, are in agreement with reference 50 if 
the species conduction currents in the undisturbed plasma are written in terms of thermal 
velocities and the above electric field and number density scales are combined to form the 
Debye shielding length for planar surfaces ( L,L,=Ap’” ). 
B. NUMERICAL DATA 

To gain further insight into the physical models presented, numerical analyses were 
conducted with data from reference 36 for 6000°K nitrogen; data for the ionization and two- 


body recombination coefficients are found in Appendix B. 


| 
(ee 
Temperature (°K), T, | 6.0000e+03 


Total Current (mA/cm?), J: used for j,;, joe 





1.0000e+02 


i 1.0000e+19 





€—————————— À 





Undisturbed Plasma Charge Density (1/m*), n,: used for n,,, A. 


Undisturbed Plasma Total Density (1/m?), N: 2.0000e+24 












Undisturbed Plasma Electric Field (V/m), E, 1.2000e+04 







Electron Diffusion Coefficient (m?/s): D/D,- n/M 3.1665e-01 


Table 4: Numerical Data for Analysis - Nitrogen 
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IV. PHYSICAL MODEL 


A. SPACE-CHARGE REGION: SHEATH 

Voltage losses at electrode boundaries, surface erosion and sheath effects play an 
important role in plasma devices and must be controlled in designs of practical interest. 
Particularly, thermal arcjets and MPD accelerators deposit between 15% and 80% of the 
input power into the anode. This presents not only a severe performance penalty, but also 
introduces thermal design problems since the heat thus generated must be radiated away from 
the thruster surfaces. [Ref. 7, 9] So a natural question: what happens at the anode? 

Near a surface (anode) in contact with a plasma there is a region, the sheath, in which 


the electric field 1s strong and in which the electron number density is much larger than that 


of the ions (see Figures (2) and (3)). 












Anode Fall 
Cathode Fall Region (see [© m 
SR Figure 3) Av, s | 
= P e: N;~D.~ Dgahe D saha 
3 £ : ; 
z S 
Y o 
EP E 
E 
Ve ES 
ONE y Sheath : Ss * Quasi-Neutral ' Plasma 
Cathode (-) Anode (+) Region 3 (Ambipolar) Region 
s Region 
E 
Figure 2: Axial Potential Profile for a High Figure 3: Qualitative Spatial 
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Figure (3) shows the qualitative behavior, in the anode fall region, of the spatial 
distribution of species densities; 1t is expected that from the wall outward, the densities 


approach one-another, eventually reaching Saha equilibrium at the plasma boundary. The 
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overall shape of the curves (Figure (3)) stem purely from the density flux = -vn-an? , 


in the absence of any other types of influences. However, it is usually the interaction of a 
systems’ many parts, not the action of each piece independently, that determines the final 
system behavior. 

The mathematical model developed for the sheath, whether anode or cathode, depends 
on the details desired: whether surface effects are considered, whether ınduced or external 
magnetic fields are included, as to how the flow chemistry is modeled, among others. For the 
cathode, sheath models usually adopt the assumptions that electrons are in equilibrium with 
the electric field and that the ionization in the sheath can be neglected [Ref. 34]. In 
comparison, the electrons in the anode fall must supply by ionization enough ions to account 
for the ion current that flows out of the positive column toward the cathode [Ref. 43: pp. 59- 
60]; electrons emerge from the column, are attracted by the anode, gain energy through the 
voltage drop and ionize neutrals. The extent of the anode sheath is then determined by the 
space charge distribution and by the magnitude of the gas ionization potential [Ref. 35: p. 
420]; hence, the associated length scale for the sheath is governed by that which results from 
the dimensional analysis of the Gauss equation (8). 

Moreover, after a current flow is established between the plates, the anode region may 
exist in a vapor that issues from the electrodes. In vacuum arcs, Miller [Ref. 44] characterizes 
the anode region as operating in one of five distinct modes, ranging from a passive, low 
current mode to a high current, fully developed spot mode. These observed instability 
phenomena are similar to those noted in MPD thrusters - anode spot formation at high current 
clearly limits anode lifetime [Ref. 6, 7, 9, 17]. What could cause such “spotting?” 

In a comprehensive review of anode-spot phenomena, Miller suggests [Ref. 44] that 
a sudden increase in the voltage noise frequently indicates that an anode foot-print could be 
forming or that a transition into the anode-spot mode could be occurring. Further, this abrupt 
change in voltage has been associated with a sudden change of ion density in the region near 


the anode, a density decrease that would leave the local electron space charge uncompensated 
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and thus produce a low-conductivity region. The sudden increase in voltage would then 
reflect this local shortage of ions - an ion depletion region or ion starvation region. Even- 
though the most suitable model for the anode region is still actively debated, the general idea 
of the appearance of an ion deficiency region near the anode as triggering the transition into 
the anode-spot mode has been fairly supported. [Ref. 44] So, to better understand the 
processes involved within the anode fall, equations (7-9) are modified as follows: 

Because the anode-spot instability is likely related to the species density distribution 
within the anode fall region, the role of net charge production, the chemistry, should be 
critical to the overall system behavior: from equations (7-9), the species densities are 
influenced by both coulombic interactions as well as by the currents and consequent 
production phenomena, the fine balance between ionization and recombination. Therefore, 
the current densities as given by equations (9a-9b) are combined with the species densities, 
equations (7a-7b): the latter are differentiated with respect to the spatial variable in this effort. 
The result gives for the space-charge region a system of two second order, ordinary, non- 
linear, coupled differential equations with one first order, coupled, ordinary differential 
equation. The non-dimensional parameters given in the bracket notation are listed in Tables 


2 and 10: 
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Essentially, the conservation laws are combined: species momentum now reflects the 
influence of charge production terms. Here, the left-hand side of equations (10a-10b) is the 
rate of change of the number density gradient, an accelerating behavior. The gradient of the 
species number density is found in the first term on the right-hand side of set (10a- 10b), a rate 
influenced by the strength of the electric field as measured by the imbalance between the 
species number densities, the Gauss equation (10c). The first term of the right-hand side of 
(10a-10b) describes the interaction of coulombic and gradient forces, the influence of the 
equation of motion as described in the earlier formulation of (1-5); the third term completes 
the formulation of motion in this model. The general concept of energy conservation begins 
to take shape with the second and fourth terms of (10a- 10b): The former has included two- 
body recombination effects as diminished not only by coulombic interactions, but also by 
diffusion. The fourth term illustrates a source of energy to the system , ionization, again as 
influenced by diffusion. Essentially, the accelerating terms on the left-hand side of (10a) and 
(10b) are described by balances of motion (gradients) and energy (ionization, recombination) 
where the exchange of energies is brought about and influenced by diffusion processes. It 
should stand to reason that as the species density gradient increases, energy losses should 
become more dominant until the gradients diminish their contribution to the flow: 
recombination should begin to take hold until eventually equilibrium is restored, balancing 
diffusion with energy source and sink. 
B. TRANSITION ZONE: QUASI-NEUTRALITY 

In order to determine the behavior of the anode fall, boundary conditions have to be 
applied to equations (7-9) or to set (10). In the case of the positive column and the cathode 
fall region, it is known [Refs. 45, 46, 47, 48] that for the isothermal case, the spatial 
derivatives of the particle density and the drift velocity become singular if the drift velocity 
attains the ambipolar sound speed [Ref. 49, pp. 77-86; 50]. This apparent singularity appears 
near the wall, marking the end of the quasi-neutral region and the beginning of the space- 
charge sheath adjacent to the wall; the ambipolar solution breaks down because the space 
charge density cannot become larger than the local density of ions or electrons. [Ref. 50; 33; 


52: pp. 202-204]. Hence, researchers have always chosen the ion sound speed as the 
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appropriate bound for the jon velocity. This condition seems a natural compromise between 
Bohm’s criterion for the collisionless sheath which requires the ion velocity exceeds or is at 
least equal to the ion sound speed and Persson’s criterion for an ambipolar flow at the plasma 
boundary, which states that the ion velocity is in fact less than the ion sound speed. However, 
these criteria are incompatible as was shown by Godyak and Sternberg. [Ref. 34] For 
example, in the numerical treatment of the plasma-wall problem without plasma-sheath 
separation, the ion velocity reaches ion sound speed ın the region where plasma neutrality is 
violated [Ref. 51: p. 51]. 

It seems that since there is no strict boundary between the plasma and the sheath, it 
is not obvious which boundary conditions work the best; much controversy remains for the 
negative electrode. Attempts have been made to bridge or to remove these apparent 
discontinuities by introducing a transition layer [Ref. 27] or by choosing various plasma- 
sheath boundary conditions [Refs. 45, 34]. For example, Godyak and Sternberg [Ref. 34] 
show that the usual choice of choosing the ion sound speed as the boundary, leads to 
discontinuities of the potential, plasma density and velocity gradients. These workers 
therefore introduced as their boundary condition a specific value of the electric field that 
provides a breakdown of the neutrality at the plasma boundary and so treated the sheath and 
ambipolar model separately, attaining a smooth transition between sheath and plasma - for the 
cathode region. In their model, however, only constant ionization was assumed for species 
production. | 

Fundamentally, the issue of transition in the cathode region remains somewhat 
controversial, a subject even more unclear when the positive electrode is considered, where 
essentially scant information is available. But a transition to ambipolar diffusion must occur 
[Ref. 52: pp. 64-68], a transition where the electric field domination shifts to that of a density 
gradient driven region, an area where now the length scale for the number density becomes 
important. 

To capture this effect and to retain the role of charge production in the dynamical 
behavior of the system, equations (7) are again differentiated allowing the species densities 


to be written in terms of the production terms. Applying the quasi-neutral condition, 
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n~n, „then gives the transition model, equation (11) which can also be obtained from the 


I e 


sheath formulation, set (10). 

















Equation (11) 


In comparison to equations (10a-10b), the space-charged model, the effects of 
diffusion are clearly seen not only on the rate of change of the number density or density 
gradient term (first term of (11)), but also on the energy balance terms, recombination 
(second term of (11)) and ionization (third term of 11). Furthermore, coulombic interactions 
are no longer evident in the energy balance, playing only a secondary role in the density 
gradient description (first term of (11)). 

A key character of equation (11) is that ambipolar diffusion is not yet fully formed, 
a consequence of the transition region in which perhaps the instability of the sheath, 
electrostatic repulsion of ions, is carried forth by gradient effects so that ions are further 
pushed from the positive plate. The growth in densities away from the plate then affects their 
free diffusion, a process evolving into ambipolar diffusion. | 
C. AMBIPOLAR REGION 

When the density of electrons and ions becomes large enough, their mutual Coulomb 
fields affect their free diffusion. This modification can be written in terms of the space charge 


field as given by: 
" Equation (6) 


From (6), it is seen that should the species flux densities not be equal, a significant charge 
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imbalance would arise. For example, if the plasma is much larger than a Debye length, it must 
be quasi-neutral and the diffusion rates of the species would adjust themselves so that the 
same rate from the plasma into the region of quasi-neutrality is achieved: the electrons have 
a larger thermal velocity as compared to the ions and tend to leave the plasma first. A 
positive charge 1s thus left behind creating an electric field of such polarity as to retard the 
loss of electrons and to accelerate the loss of ions: the plasma, by its very definition, is a 
neutral fluid. [Ref. 38: p. 139] From this formulation, the ambipolar diffusion coefficient can 


 D,u* Dy, 
W+, 


be determined as D 


a 


where the Einstein relation may be applied to write 


mobilities in terms of diffusions. It should be pointed out that the ratio D/u represents a 
measure of energy [Ref. 52: pp. 60-63], and that if T,-T; , as 1s assumed in this work, the 
ambipolar electric field tends to enhance the diffusion of ions by a factor of two, a diffusion 
that is primarily controlled by the slower species: D,~2D,. The primary function of the 
ambipolar field is then to serve as an energy-exchange mechanism for transferring random 
kinetic energy from the electrons to the ions such that the ion diffusion continues away from 
the positive plate. [Ref. 45] 

The model for the ambipolar region is thus formulated, again following naturally from 
equation (10) or from equations (7-9), the sheath model: proceeding as before, equations (7) 
are differentiated with respect to the spatial variable allowing set (9) to contribute to the 
densities. However, unlike previous, now set (7) is multiplied by the respective species 
diffusion coefficient, the result which gives one second order, ordinary differential equation 
for the density in terms of not only the necessary chemistry, but also in terms of the ambipolar 
diffusion, D.; 
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nv Equation (12) 


In comparison with equation (11), it is immediately seen that the density gradient 
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contribution, the rate of change of number density, is no longer evident. The only terms 
remaining in (12) are those of charge production; the system chemistry. As the flow 
approaches equilibrium near the plasma boundary, the number density must approach that of 
Saha so that the definition of “plasma” is upheld. 

The physical model in summary: the electric field pushes the flow from the electrode 
surface toward the density gradient dominated region in which the density gradient pushes the 
system toward the plasma boundary where eventually energy conservation must be satisfied. 
Ionization must balance recombination such that Saha equilibrium is attained at the plasma 
boundary. It is in this region that the largest of the length scales has influence: the current 
or ionization scale. This last scale defines the extent of the anode fall region, a region where 
flow chemistry becomes paramount, determining plasma equilibrium. 

But is the anode fall stable? Previous attempts [Ref. 53] at numerical analysis have 
proven challenging at best. So a better understanding of the nature of equations (7-9) and 
sets (10-12) is needed: in the next chapter, non-linear analysis techniques are applied to 


mathematical models presented previously. 
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V. ANALYSIS 


A. LOCAL ANALYSIS 
To begin the analysis of the anode fall, non-linear methods as briefly outlined in 
Appendix C are used: further information on particular ideas can be found in references 54 
and 55, among others. 
jr Space-Charge Region: Sheath 
To begin an analysis of either equations (7-9) or equations (10), non-linear methods 
as outlined in Appendix C are used: first, equations (10) are written as a system of first order 
differential equations, then stagnation points within both systems, equations (7-9) and the first 
order formulation of equations (10), are determined. Next, the local behavior about these 
stagnation regions in terms of system characteristics is evaluated, leading to a qualitative 
stability picture of the system dynamics. 
a. Equilibrium (Stagnation) Regions 
First equations (7-9) and equations (10), the anode sheath models are 


analyzed (non-dimensional parameters are given in Table 2): 





























KR) [ 1-2 1,17: d(ñ,) -z [ul 
I = p- en a = +id > 1 ds 
dy El ni b| ay Ic |n E | | D. Equations (7a)-(7b) 
+ -[e |z, -|y D Equations (8) 
d(j,) ENT Gia es 
d? - -|e va, * |^ I E E va, - | p laa, Equations (9a)-(9b) 
where D, ‚=D, NM 
u (en. ee — 2 
dy = a E T + D CH A + lae n; - ne Equation (10a) 
d*(n,) r d(n,) dh —— — dg i=- 
ay E [el dv + Z -ce i + le I = nio Equation (10b) 
> x l n - ls. Equation (10c) 


25 


The original model, equations (7-9) can immediately analyzed, whereas equations (10) 


have to be recast as a set of coupled, first-order differential equations. 








i a ae jr - ED i 
dy D, n E Equation (10a) 











Equation (10b) 


lus vi. Equation (10c) 


To find the stagnation points within these first-order system, the derivative terms are each set 
to zero and non-dimensional parameters from Table 2 are applied. The results give for the 


equilibrium criteria: 
1. From equations (8) and (10c), n," -n,'-n' as expected; this result is then 


incorporated into equations (9a-9b) as well as (10a-10b). 
2. Similarly, from equations (9a-9b) and (10a-10b), n*=0 or n*=—“— 3a 


stagnation point appears when n,'=n"=n"=0 and when the plasma boundary 


is approached. The latter condition immediately emphasizes the flow chemistry, 
ionization and recombination. Without recombination, the second stagnation would 
not be present, as seen from either equations (9a-9b) or equations (10a-10b) with 
d/dy[ ]-0. The former equilibrium, on the other hand, seems artificial at first, unless 


the species densities behave such that they cross when n,'-n'/-n'-0 , 
immediately implying a quasi-neutral condition since equation (8) or equation (10c) 


se 


require 7,'-n '-n' at stagnation. In general, the fixed points found stem purely 


from the chemistry behavior, equations (9a-9b) and consequently (10a- 10b) in which 
ionization and two-body recombination are the only processes included. If, on the 
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other hand, three-body recombination is added to the equations, then additional 
stagnation regions become evident since equations (9a-9b) are now cubic. The 


equilibrium region near 7,’ =n," =n*=0 still remains, however. Similarly, if 


ionization is the only production term considered in the sheath model, then the 
stagnation region near the plasma boundary “disappears.” The condition 


n'-n'-n'-O0 remains for all cases: ionization only, ionization with two-body 


recombination, ionization with two- and three-body recombination. So it appears that 
continuity (equations (9)) requires flow equilibrium when perhaps quasi-neutral 
conditions are initially met, at the sheath boundary, and when the plasma interface is 
approached. It is recombination that establishes the latter while the presence of 
ionization stipulates the former. 


3. The current densities at the equilibrium position can be found by using n*=0 
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since n,=n,=n, when n,=n,=n" . Here, the influence of not only 


e 


ionization and recombination, but also of diffusion on the current density distributions 


are clearly seen. Similarly, from equations (7a-7b), ™ SOM VD 


4. Charge conservation is used to isolate the value of the electric field at stagnation 
vv 


near the plasma boundary, where  m'- That is, 





en, 


VJ-20 —J-j *j,-J or DE gives 
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from which E is obtained graphically since the right-hand side of this expression 
simply represents a constant. Applying the Arrhenius function [Ref. 52: pp. 94-109] 


ZA 


E V 


=— and using numerical data for nitrogen 


V 
o 








for the ionization, v-e 


yields E'*«0.54 .Fromthis, —j, «0.99816, j,' «0.0005445, n «0.15695 


These results show the relative magnitude of the electron and ion currents within the 


anode fall, j.»j; as expected. If, on the other hand, Tre '20 Vn'-0 isused, an 


inconclusive result is obtained: E” cannot be isolated from either equations (7-9) 





-1 
or equations (10). Further, v-elf | —. does not permit zero field. Therefore, 
NS 


it seems reasonable that a finite field should exist over the entire range of the anode 
sheath. Von Engel [Ref. 56: pp. 12-20] highlights such arguments: a finite field must 
be allowed, which could then introduce a potential trough. A virtual anode is thought 
to exist 1f the space charges have the same sign as the nearby electrode, with the 
emitting surface located at the potential minimum [Ref. 56: pp. 14-17]. 


The fixed points found in the model represented by equations (7-9) and in part by set (10) are 


listed in Table 5 below, featuring the role of ionization, recombination and diffusion in the 
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existence of the fixed point near the plasma boundary. Again, continuity (equations (9)) 
seems to require flow equilibnum when perhaps quasi-neutral conditions are initially met, the 
sheath boundary, and when the plasma interface is reached. It is the recombination 
phenomena that establishes the latter stagnation region while the presence of ionization 


stipulates the former. 


T 


n, 


Near the Plasma Interface Numerical 
Value for 


Nitrogen 


0.15695 


0.54 
i GRE: gEV D, 
KT ,0j, KT Oa), 
0.99816 
0.0005445 


Table 5: Stagnation Regions in General Sheath Model - Ionization with Two-body 
Recombination 
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b. Jacobian at the Equilibrium Positions 
To determine the system characteristics about its fixed points, a local Taylor 
series expansion is carried out about the equilibrium of interest (Appendix C); the rate of 


change of the flow is given in terms of the Jacobian, which for equations (7-9) becomes: 














i ; Jo] 
[a] E 0 [a]7, 0 5 
0 -[c]E -[e]n, E 0 
Å inner = [e] -W 0 0 0 
(ala, [given 2 0 
E 
E = A [i]r,v 
Em [¿Jv-[7]n, = 0 0 


Here, six variables affect the system dynamics: ion and electron number densities, 
electric field, species diffusivities, and ionization. Of these, the diffusion coefficients are 


assumed constants, so essentially four variables drive the behavior: species densities, electric 


field and consequently ionization. If Ai is evaluated at the stagnation points, the 


resulting eigenvalues describe the dynamics in the neighborhood of the fixed point (Appendix 
C). Specifically, at the plasma boundary equilibrium with nitrogen data: 
A, = 1.376704 4, =0 


A, =0.630455 A,=-1.37437 
A. = -0.63279 


For this case, the eigenvalues are real and distinct, leading to a multi-dimensional type of 
“saddle” behavior (Appendix C). The stagnation region near the plasma boundary 1s unstable: 
the largest eigenvalue is positive, indicative of exponential growth. Appendix D lists the roots 
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ofthe characteristic polynomial of A. . It appears A, gives exponential growth - the 


inner 


process of ion repulsion from the positive plate forces the instability of the sheath. The sheath 


drives the system toward the first stagnation criterion, toward quasi-neutrality. 


Provided a finite field exists when n,'=n"=n"=0 „then the fixed point in its 


vicinity gives a similar qualitative behavior: the eigenvalues of the resulting Jacobian are given 


by 





a 
Age [c]E , E where v~e z|. v 
; 2 2 V 


o 





These eigenvalues are always real with two characteristics equal while the remaining are 


distinct: Table 6 lists some numerical examples for nitrogen data. 


D 


10.00 


-2.9986 | -9.99947 
J -0.0014 | -0.00053 


Table 6: Eigenvalues of Stagnation Region n=n,=n=0 with Finite Field - Numerical Values 
for Diesen 





Again, the essential feature of the sheath is instability: ions are driven from the 


positive wall as required by electrostatics. The stagnation region near n'*=n*=n*=0 is 
unstable, perhaps pushing the system away toward the recombination induced equilibrium 


near the plasma boundary. It is the strength of the initially applied electric field that imposes 


the stability characteristics for the fixed point n,*=n*=n*=0 . This leads to the 


following question: is the sheath model, as presented by either equations (7-9) or set (10), 
an example of a system in which a certain type of “fixed point” is always present, a stagnation 
region satisfying quasi-neutrality and thereby representing the sheath boundary? If so, does 
this equilibrium present a type of transcritical bifurcation (see Appendix C) whose stability 
perhaps depends on the initially applied electric field? The foregoing does suggest so. 

The precarious nature of the fixed points found in the sheath model suggest that 


numerical integration methods will not prove fruitful over the entire domain of the anode fall. 
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However, by exploring the various regimes of equations (7-9) in the form of the second order 
differential equations given by equations (10), further insight into the anode fall could be 
expected. 

2 Transition Zone: Quasi-Neutrality 

Previously, it has been shown that the sheath is inherently unstable, that the ions are 
electrostatically repelled from the electrode (anode) surface. Ifthrough that process quasi- 
neutrality occurs, then for steady-state conditions either density gradients or diffusive 


processes have to drive the system toward the plasma boundary since the condition of 


nn, stipulates the electric field be constant; so something else has to continue driving 


l 


the ions away from the anode wall. Equation (11) highlights the density gradient and is 


considered next as a set of coupled first-order differential equations: 
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Equation (11) 


where the electric field is constant due to quasi-neutral conditions, n-n, 
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a. Equilibrium (Stagnation) Regions 
To find the stagnation points within these first-order system, the derivative 
terms are each set to zero and non-dimensional parameters from Table 2 are applied. The 


results gives for the equilibrium criteria: 


1. Setting d [ ]/dy-0 in equation (11) immediately gives the previously found 


vV 
2 The fixed 





stagnation condition near the plasma boundary: n"= 


o no 


point n*=0 again seems to present a mathematical phenomenon since equation 
(11) does not describe the sheath condition. If quasi-neutral arguments are again 
used, then perhaps the condition m*=0 may have meaning - the plasma shields 
itself against external influences (electric field), and so, quasi-neutral characteristics 
are fundamental to its definition. The stagnation regionat n*=0 does, perhaps, 


present a physical phenomenon: the onset of quasi-neutrality or the sheath boundary, 
found as a consequence of ionization present in the sheath model and in equations 
(11), the “transition” model. 


2. The magnitude of the constant electric field, and therefore constant ionization, at 
the start of the quasi-neutral transition is determined by the matching conditions when 


n~n, . These criteria are either found through analytically matching equations (10) 


3 e 


with equations (11) for some small parameter(s), or numerically. Numeric methods 
will be considered later in this work. Suggestions toward an analytic matching 
attempt are provided in paragraph B below. 


b. Jacobian at the Equilibrium Positions 
To determine the system characteristics about its fixed points, a local Taylor 
series expansion is carried out about the equilibrium of interest (Appendix C); the rate of 


change of the flow is given in terms of the Jacobian, which for equation (11) becomes: 
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Since the electric field becomes constant at the onset of quasi-neutrality at the sheath 


boundary, the electric field thus locks in the behavior of the transition region. To determine 


these stability characteristics, the general eigenvalues of A are determined: 


transition 





; _aE(D,-D,) a?E’(D,-D,) -4¥(D,+D,)(bi+dg) +8n(D,+D,)(bj + dh) 
t2  2(D,+D,) 2(D, * D,) 


where the non-dimensional parameters are listed in Tables 2 and 10. These eigenvalues 


determine the dynamical behavior of the transition when Á is evaluated at a 


transition 


fixed points, albeit the electric field strength is determined from the sheath as the system 


reaches 7-n, . By inspection of the first term found in these eigenvalues, the following 


is observed: 


1. Ifthe field is positive, then the system displays exponential decay, stability; but 
electrostatics demands that ions be repelled from the positive plate. This leads to the 
expectation that the transition zone is also unstable, similar to the sheath. 
Consequently, the field should be negative at the point of quasi-neutrality. 


2. Using the stagnation criteria for n*=0 with finite field indicates that the 


characteristics should display stable behavior for positive field, the diffusion difference 
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(D,-D,) is negative, and unstable behavior for negative field. The initially applied 
electric field strength then does determine the stability of the fixed point n*=0 


? 


an equilibrium which appears when ionization is present in the model. 


Vo V 





3. To find the behavior near the plasma boundary when  n'- , nitrogen 


œ, n 


20 


data along with an arbitrary selected field strength, E (;_; ) = 5.00 are used to 


give the eigenvalues A,=-0.041 and A,=5.026: the eigenvalues are real, distinct and 
predict a saddle-type behavior near the plasma boundary - a robust linearization (see 
Appendix C). For this case, the flow field, locally, is characterized by a very strong 
unstable manifold corresponding to A,. A similar result was found by analyzing the 
sheath model presented by either equations (7-9) or equations (10). 


4. On the other hand, if the sheath behavior at quasi-neutral formation gives a 
positive field, e.g., E (y..- ) =5.00 , then the Jacobian eigenvalues predict stability 


locally about the fixed point: A,=0.021 and 1,=-5.013 where now the saddle behavior 
is marked with a strong, stable manifold, contrary to the sheath analysis earlier. 


Thus, the criteria found at the onset of quasi-neutrality at the sheath boundary form 


the initial conditions imposed upon the transition zone and hence govem the stability 


vv 
o 





characteristics of the fixed points at n*=0 and n*= 


A, AN 


2 0 


3. Ambipolar Region 

In the sheath, ions are repelled electrostatically until quasi-neutrality occurs, where 
chemical processes balance such that the ion density approaches that of the electrons - the 
electric field becomes constant. If this balance does not reflect total charge neutrality, i.e., 
Saha equilibrium, then other processes (the field is now constant) must drive the system 


toward such equilibrium; however, a stabilizing mechanism must be available once the 
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densities are near the plasma boundary. So as the species density gradient increases, energy 
losses should become more dominant until the gradients diminish their contribution to the 
flow: recombination should begin to take hold until eventually equilibrium is restored, 
balancing diffusion with energy source and sink - ambipolar diffusion becomes the vehicle 


which takes the flow density toward the plasma stagnation. To translate these physical 


processes to a mathematical model, quasi-neutral conditions, n,-n, ‚are applied to the 


sheath equations (10) from which equation (12) results naturally, highlighting the ambipolar 
diffusion in its full form: 





dam (61. dhl= [bias] e 
^ & 2 D, D, nv Equation (12) 
or 
, Un) _ x, 
dy 
dx) _ Y a A EXE Equation (12) 
dy D, D, D, D, 


where the electric field is constant due to quasi-neutral conditions, n,~n, . 


a. Equilibrium (Stagnation) Regions 
To find the stagnation points within this first-order system, the derivative 
terms are each set to zero and non-dimensional parameters from Table 2 are applied. The 


results gives for the equilibrium criteria: 


1. Setting d [ ]/dy=0 in (12) immediately gives the previously found stagnation 


* O 


condition near the plasma boundary: n* = 
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2. Similar to the transition zone, the fixed point n*=0 again presents a 
mathematical phenomenon until physical arguments are applied as before. 


n*=0 perhaps indicates the onset of quasi-neutral conditions. 


The magnitude of the constant electric field, and therefore constant ionization, at the 
start of the quasi-neutral ambipolar region continues from the transition zone and again locks 
in the behavior of the ambipolar region: the initial field strength determines the initial 
conditions for, and consequently the stability of, the transition zone. The way in which the 
density gradient then develops influences the initial conditions to the ambipolar region. All 
in all, the initially applied electric field strength impacts the stability of the anode fall, in 
entirety - the electric field directly influences ionization and hence the stagnation found near 
the onset of quasi-neutrality. 

b. Jacobian at the Equilibrium - Analysis 
To determine the system characteristics about its fixed points, a local Taylor 
series expansion is carried out about the equilibrium of interest (Appendix C); the rate of 


change of the flow is given in terms of the Jacobian, which for equation (12) becomes: 
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with general eigenvalues 





te 


If the stagnation conditions near the plasma boundary are applied, it can be easily be seen that 


the ambipolar eigenvalues are always real and distinct, again predicting the saddle-type 
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behavior (Appendix C). In comparison to both the sheath and the transition zone, the electric 
field no longer is found explicitly part of the Jacobian eigenvalues, however an implicit role 
is still evident in terms of the ionization. Furthermore, diffusion and chemistry are part in all 
eigenvalues, from sheath to ambipolar. 


Proceeding as before, the Jacobian is evaluated at the plasma boundary stagnation, 


«| 


ek y 
n = 
e n, 





, using nitrogen data. The initial field strength is chosen arbitrarily; “initial” is 


defined as that condition found at the onset of the ambipolar region as determined by 
matching criteria between sheath and transition, at the beginning of quasi-neutrality. The 
resulting Jacobian eigenvalues or characteristics are again robust (system non-linearities do 


not change the qualitative nature of the flow ) [Ref. 54: p. 163]. That is, if 


E(y-_-)=-5.00 ,À,57£ 37.13, whereas if. E(y. - ) 35.00 ,4,,72530.4. Although 


the eigenvalues are real and distinct, they are of the same magnitude: exponential growth is 


exactly offset by exponential decay - a hallmark of energy conservation or homoclinic 


behavior (see Appendix C). The phase portrait for equations (12), — n(y) vs dp à 
J 


using nitrogen data, is shown in Figure (4). 


39 





d[n(y)]/dy 










) f P 4 
ty LIN 


E 
Quasi-nentrality: ‘condition 
Pede O bit) f 





WW / 


Region [ 






n(y) 


N 
N 


Figure 4: Phase Portrait for Ambipolar Region, Equation (12), with Nitrogen Data for 64 

Different Initial Conditions using Ionization and Two-Body Recombination 

Can the stagnation near the plasma boundary be reached? Is this region stable? The 
behavior of the ambipolar region is driven by the initial conditions as applied to equations (12) 
as shown by figure (4): if the initial conditions are such that region I is dominant, then the 
manifolds of the stagnation near the plasma boundary, the saddle behavior (see Appendix C), 
confine the system trajectories. If, on the other hand, the initial conditions are such that 
region II is dominant, then the system will become unstable. The manifolds of this 
linearization form the threshold, at least locally, which perhaps determines the requirements 


for quasi-neutrality - the homoclinic trajectories, as given below: 
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Here, the contributions of ambipolar diffusion, ionization and two-body recombination are 
clearly seen; two-body recombination is found in the non-dimensional parameters [h] and [j]. 
More often than not, these homoclinic trajectories can be volatile: any numeric round-off 
error incurred through whatever computational algorithm used, will induce sufficient non- 
linearities that will prevent the solution from being “pinned” to these manifolds, the manifolds 
will break. Even if it were possible to achieve this exact trajectory, the initial conditions 
imposed by matching the transition region to the ambipolar equation will most likely make 
such efforts challenging. That is, the initial conditions provided to the quasi-neutral region 
by the sheath determine the behavior of the ambipolar evolution: either stable periodicity as 
shown by region I in Figure (4), or exponential growth as indicated by region II. In either 
case, the behavior of the sheath ultimately determines the behavior of the outer regions and 
therefore the anode fall, in general. 
4. Comments/Observations 


The preceding suggests that: 


1. Ions, created through electron bombardment of neutral atoms throughout the 
anode fall region, are repelled electrostatically from the positive electrode toward a 
location of quasi-neutrality at the sheath boundary as well as toward an unstable 
stagnation region near the plasma interface. These equilibria form as a result of the 
chemistry present. That is, continuity (equations (9)) seems to require flow 
equilibrium when quasi-neutral conditions are initially met within the sheath, and when 
the plasma boundary is reached. It is recombination that establishes the latter while 
the presence of ionization stipulates the former. Fundamentally, the sheath as 
formulated is inherently unstable - ions are driven out of the system, toward the 


4] 


negative electrode: it is the role of the anode, after all, to provide the positive current 
- to the cathode [Ref. 43: p. 59]. 


2. At the moment of quasi-neutrality, when 75-7, , there appears a conservative 


phenomenon, such that the energy growth is matched by the energy loss - a chemical 
balance. Is the sheath model as presented, an example of a system in which a certain 
type of “fixed point” is always present, a stagnation region satisfying quasi-neutrality? 
If so, does this equilibrium present a type of transcritical bifurcation whose stability 
depends on the initially applied electric field? The brief analysis given in this chapter 
does suggest such behavior - the onset of quasi-neutral characteristics, the sheath 
boundary, might be indicative of a transcritical bifurcation. 


3. The quasi-neutral balance, however, is insufficient to provide the necessary 
stability condition to ensure the system progresses toward, and remains in the close 
vicinity of, the plasma boundary fixed point. As such, any numeric solutions will most 
likely be difficult to achieve. That is, the initial conditions provided to the quasi- 
neutral region by the sheath determine the behavior of the ambipolar evolution: either 
stable periodicity as shown by region I in Figure (4), or exponential growth as 
indicated by region II. In either case, the behavior of the sheath ultimately determines 
the behavior of the outer regions and therefore the anode fall, in general. 


B. ANALYTICAL MATCHING - SUGGESTIONS 

The analytic analysis presented in this chapter suggests that the anode fall model used 
in this work does not represent sufficient dissipation, preventing the system from adjusting 
itself, to adhere, to the “stable manifold" of the stagnation region near the plasma boundary. 
The model is most likely incomplete in its physical representation. Perhaps analytically 
matching the second order sheath equations to the transition set which are then matched to 
the ambipolar equation, may reveal an inconsistency between some of the physically 
represented terms. Albeit this process seems challenging at first, perhaps some of the ideas 


suggested next could be applied: 


l. Three length scales have been found when equations (7-9) were non- 
dimensionalized (Appendix A and Table 3). Three “small” parameters could then be 
used to match equations (10-12) as follows: Assume the sheath equations (10) 
behave rapidly over a very small distance, hence a fast changing parameter should be 
applied, 1/6=O(L,). Further, assume the ambipolar model (12) as varying slowly in 
comparison, €~O(L, or L;). In between is a region of transition where n=e/0. But 
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where are these equations matched? 


3. To determine the location of matching, the fixed points found from equations (7-9) 
could be used as a starting point. The linearization of equations (10-12) about these 
stagnation regions, then provides the criteria for matching: the respective manifolds 
describe the number density and the density gradients evident near these fixed points. 
The interpretation: the manifolds given by the linearization about the fixed points 
gives the conditions for matching, criteria which must be met by the preceding region 
as the system is driven into the next zone. The manifolds for the fixed point near the 
plasma boundary have already been presented in section A.3. above. The manifolds 


for the fixed point near the sheath boundary can be found from A RE DY 


using any symbolic mathematic utility such as MAPLE or MATHCAD. The 
following is an illustration of the above concept: 
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3. The initial conditions that are used for the sheath, however, remain a “shooting’ 
problem. Perhaps a numeric iteration scheme where many initial conditions are tested 
might reveal some pattern for better selecting these conditions. 
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VI. NUMERICAL RESULTS FOR NITROGEN 


Using available MATLAB 4-5-stage Runge-Kutta algorithms, both equations (7-9) 
and equations (10), as first order differential equations, are integrated numerically. The data 
thus obtained were then exported into a Quattro-Pro spreadsheet for plotting and easier 


analysis. 
1. Length scales (Table 3) to determine the non-dimensional parameters: 
kT Jis 
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2. Non-dimensionalized initial conditions (Table 4 lists the numerical data for 
nitrogen): 
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The dnft currents are used to generate the appropriate initial conditions at the 
electrode wall since at the wall, the field induced motion will exceed that induced 
thermally; it is assumed that initially any density gradients are zero. To gain insight 
as to the impact of initial conditions on the system, the initial electric field 1s varied 
from 0.01 to 20(10*) V/m. The results are given below. Of note, the method used 
to non-dimensionalize the variables (see Appendix A) gives rather small “non- 
dimensional numbers” in terms of magnitude, especially when the integration results 
shown are “close” to the electrode wall. 


45 


A. SPACE-CHARGE REGION, SHEATH 


In overview: For ŒE (0)< E, »reshoig ^ 0.04 „the ion density behaves as expected, 


increasing monotonically, being primarily influenced by its gradient. Both the low initial field 


strength and the ion density gradient are insufficient, however, to push the species toward 


quasi-neutrality when nn, When E(0)2E, ,,,-0.04 , the initial field becomes 


strong enough to affect the ion density as compared to its gradient; the ions are pushed 


toward what appears as a density trough. Furthermore, for certain initial values of 


An=|n,-n,| vs E(0)  , the ions are pushed past their minimum density point, 


eventually approaching the electron density distribution, after which both densities remain 
close in magnitude for a short distance; essentially, for adequately strong initial electric field, 
quasi-neutrality is formed for a brief distance. The quasi-neutral characteristic 1s inherent 
within the energy balance found in the ions, the balance between ionization and 
recombination; without this balance, the space-charge region remains. 

1. Varying Initial Conditions: Impact on System Behavior 

To gauge the impact of various initial conditions, specifically combinations of species 
density and electric field, on the sheath dynamics, equations (10) were integrated for 
n,(0)=10* [1/m?], n.(0)=10* [1/m?] and E(0)=6(10*) [V/m] using the length scales as given 
in Table 3. These conditions were chosen after varying the initial conditions as will be 
described below. 


Specifically for the ion density in the sheath: 
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lon Density vs Distance from Wall 
Varying Initial Electric Field 
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Figure 5: Ion Density vs Distance - E(0) < E mreshotd 


Figure (5) shows the ion density vs distance from the electrode wall for distances very 
close to the electrode surface - in comparison, see Figure (8) given later in this chapter. 
Again, the non-dimensionalization used makes the magnitudes for the density rather small, 
but these magnitudes represent non-dimensional quantities and so should be multiplied by n, 
(see Table 4). 

At the electrode wall, in its immediate vicinity, there are at least two major vehicles 
which drive the system away from the wall in steady-state: the initial electric field strength and 
the density gradient, such that if the field strength is low enough, the gradient dominates the 
system whereas if the field strength is high enough, the gradient influence becomes secondary 
in nature. Figures (6)-(7) depict the magnitude of the ion density gradient vs distance from 


the electrode wall and the electric field vs distance, respectively. Comparison with Figure (5) 


leads to the speculation that for low electric field strength, E(0)< E, ee 
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ion density is governed by its gradient, whereas when E (0)> Enresnoia=0-04 ‚the electric 


field drives the ions from the wall. When E (0)=0.04 the boundary is reached where 


initial density gradient effects are surpassed by those of the field (Figures (5) and (6)) . As 
shown in Figure (5), the ion density appears to increase from the electrode wall for low initial 
field; yet this growth does not continue for very long. Essentially, the field strength is 
insufficient to drive the density further while the density gradient is too weak to provide the 
necessary impetus for the ion and electron densities to match. In general, Figures (5)-(7) 
show results for integrated variables in close vicinity of the electrode surface - in comparison, 
see Figure (8). 
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Figure 7: Electric Field vs Distance: E(0) < Egreshold 
A density minimum occurs whith the initial electric field strength increases (Figure 
(8)). Moreover, the density decline reaches a minimum according to the electric field strength 
imposed. A constriction is suggested where the location and magnitude of these minima 


appear affected by the initial field strength as shown in Figure (8). 
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Figure (8) suggests the extent of the sheath as order 10($), in agreement with the 
theory which requires the shielding to be of Debye length: Ap-[e,kT,/q'n,]^. But only 
sufficiently strong initial electric field intensities are able to push the ion densities past their 
wells, past their constriction (Figure (8)). 

In comparison, the electrons in the anode sheath behave as expected, being less 
influenced by electrostatic repulsion than the ions. Figure (9) shows that for very low initial 


field strength, the electron density increases monotonically, but only very slightly. 


Electron Density vs Distance from Wall 
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Figure 9: Electron Density vs Distance - E(0) < Enreshold 


So, it is expected that increasing initial field strength should give rise to increasing 
electron density - this is the general case (Figure (10)). The electron density does decrease 
very slightly before the location of ion minimum is reached as illustrated by Figures (10)-(11); 
Figure (11) represents a more detailed view of Figure (10) in the region before ion minimum. 
Similar to the ion density, only sufficiently large initial field magnitudes push the electron 
density past the location of ion well (Figure (10)). From a physical standpoint, it should be 
remembered that both species are affected by electrostatic forces: the ions, however, are 
influenced to a much greater extent due to the electrostatic repulsion required when like 


charges are sufficiently close; qualitatively, these are the observations found in the preceding 
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analysis (Chapter V). 
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If the initial field is sufficiently strong to push the ion density past its minimum, the minimum 


occurs when the electric field curvature changes: that is, the field curvature as written in terms 


"i 
of its potential (7) changes signs. Shown in Figure (12) are the results for 
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|, n,(0)=10% T 


E(0)=5.08 with n,(0) = 10° A 
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Figure 12: Densities & Electric Field vs Distance - E(0) > E, 4,475.08 


Although a quasi-neutral region forms when the initial electric field is sufficiently strong so 


as to push the ions past their well, the number density near the plasma boundary is not 


reached: n*=0.15695 as compared to the maximum density achieved 


when n,-n, where n isontheorder of10°. Since n,~n, requires E = constant, 


something else is needed to drive the system toward the plasma stagnation - some form of 
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transition. To summarize the effect of initial conditions: 


It appears that the anode fall is comprised of diverse behavior, in part determined by 


the initial conditions chosen: there appear ranges of An- |n; -n,| vs E(0) for which 


both equations (7-9) and equations (10) naturally form a region of quasi-neutrality whereas 


forotherranges of An=|n,-n,| vs E(0) ,there seems insufficient initial energy to 


push the ion density past its minimum. To gauge this effect, equations (10) were integrated 
again but now E(0) was varied from 0.01 to 20(10*) [V/m] for fixed An ranging between 10? 
to 10* [1/n*]. Figure (13) shows the results in semi-logarithmic format: the data points are 
plotted for which E(0) is sufficient to drive the ion density past the minimum, toward the 
electron density. From this plot, a type of cut-off" region seems to form: to the left, the ions 


remain locked in their density well, to the right, the ions are driven past this minimum toward 


n~n, 
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For the region of initial conditions for which n,-n, is reached, an interesting 


result appears when the ion density is plotted against the corresponding ion density: a type 
of phase portrait for the ions. Specifically, a multi-dimensional homoclinic-behavior results 
(see Appendix C), a characteristic commonly associated with energy conservation: for 
sufficient initial field, energy requires a balance between gains and losses, between ionization 


and recombination. So, if the initial conditions are such that the ion density is pushed past its 


well toward quasi-neutrality, the event of n,~n, closes the homoclinic orbit and energy 


is balanced: a comparison between the case E (0)=5.00 and E (0)=4.17 is shown in 






Figure (14). 
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Figure 14: Anode Fall with Sufficient Initial Field - Phase Portrait for Ions 


The integrated results of equations (10) not only provide the species density behavior, 
but also their gradients. As such, a natural question arises as to the behavior of species drift 
and random currents: the behavior of the species under the influences of an electrostatic field 


(potential) as compared to the behavior in an environment free of potential. 
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2: Drift and Random Currents 
The species densities are not constant throughout a region which is influenced by an 


electric field; a uniform electric field applied to a cloud of constant charge q will produce a 


directed motion along with a variation of the concentration. From this principle, the species 


drift currents are found [Ref. 57: pp. 53-55]: 
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where the Einstein relation for mobility has been applied. To obtain the random current, the 


thermal velocity 1s used such that 


— 1 PEN 1 = 2k T . 
| 1 Qn, Pie Vtnermaı,, | 1 mud E Equation (14) 
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The effect of varying initial conditions, An vs E(0), aré again determined in similar procedures 


as for the species densities above as shown next. First, the 10n currents are considered: 


, the random ion current increases, 


For field strengths E(0)<E, ,,=0.04 


following the tendencies of the ion density (see Figure (5)). When the initial field is increased 


, the random and drift currents intersect as a function 


such that E (0)> E,, eshold = 0-04 


of the rate of decrease or strength of initial field. Again, the random current (equation (14)) 
follows the density distribution for isothermal conditions. Figure (15) summarizes these 


results for locations “close” to the electrode surface. 
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Figure 15: Random & Drift Ion Current Density vs Distance - E(0) < Ennreshod 


If the initial electric field 1s sufficiently large so as to push the ions past their density 
well, the drift current follows accordingly. Furthermore, the random current follows the ion 
density as shown in Figure (16). 
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Of note is.that ion both the drift and random currents. seem to approach-a-limiting 
value-as the realm. of quasi-neutrality is approached; perhaps the saturation current as 
theorized by Child and Langmuir [Ref 58: p. 238]. Next, the electron currents are 
considered: 

For sufficiently large initial field, the rapid variation in the electron density gradient 
(see Figure (11)) is markedly observed in the corresponding electron drift current - Figure 
(17). Additionally, the electron random current does not intersect the electron drift current, 
in comparison to that of the ions. Instead, the random electron current follows the tendencies 
of the electron density (Figure (18)) where the negative charge of the electron reverses the 
magnitude of the electron density plot (see Figure (11)), in comparison to the ion results. 
Furthermore, when a sufficiently large initial field is applied, the electron drift current 
increases past the point (spatial coordinate) of the density minimum, diverging spatially from 
the random current (Figure (19)). 
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B. NUMERICALLY MATCHING 
Although a quasi-neutral region forms when the initial electric field is sufficiently 


strong so as to push the ions past their well, the number density near the plasma boundary is 


not reached: n*=0.15695 as compared to the maximum density achieved 


when n,-n, where n  isonthe order of 10”. Since n,~n, requires E =constant, 


something else is needed to drive the system toward the plasma stagnation - some form of 


transition, for without it, the species densities diverge (Figure (20, 21)). 
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Figure 20: Anode Fall - Without Numeric Matching 
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Figure 21: Anode Fall - Without Numeric Matching (Expanded View) 


In steady-state with isothermal conditions, with the electric field constant when 


n,-n, _ , there remain two primary phenomena which could drive the density toward the 


plasma equilibrium: density gradient or ambipolar diffusion. To obtain the overall effect, first 


the sheath equations as given by set (10) are matched numerically to the final ambipolar 


equations, (12) using for the matching conditions n,-n, that result from equations (10) 
provided the initial field is sufficiently strong to allow the ions to proceed past their density 


well; the results for E(0)=6.1(10*) [V/m] and n;(0)=10° [1/m?] with n,(0)=10'? [1/m?] are 
shown in Figure (22). 
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Figure 22: Anode Fall - Numeric Matching without Transition 
The density remains below the plasma boundary stagnation (Figures (22, 23))- 
ambipolar diffusion is insufficient to drive the density toward the plasma boundary. The 
cosinusoidal behavior shown is determined by the matching conditions imposed by the sheath 
equations (10) upon the ambipolar equation (12) - Chapter V. 
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Figure 23: Anode Fall - Numeric Matching without Transition (Expanded View) 
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If, on the other hand, a full transition region is considered in which the density 
gradient is the primary source to drive the system, then the plasma boundary density is 
approached. The transition equations do not, however, reflect the impact of the density 


growth on free diffusion - the formation of ambipolar diffusion. Hence, using the plasma 


boundary stagnation, n*=0.15695 „along with the electric field intensity at the onset of 


n.-n, as initial conditions (or matching conditions) to equations (12), the transition 


equations (11) are numerically matched to the ambipolar formulation, set (12). Figure (24, 


25) illustrates the results for E(0)-6. 1(10*) [V/m] and n(0)-10 [1/n? ] with n,(0)=10'° (1/nT]. 
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Figure 25: Anode Fall - Numeric Matching with Transition (Expanded View) 


The ambipolar region shown in Figures (24, 25) clearly shows the very sensitive 
nature of this zone; the initial conditions to the transition region the consequent matching to 
the ambipolar zone determine the characteristics of this latter region (see Chapter V, Figure 
(4)). Specifically, the behavior of the anode fall region, as modeled by equations (7-9) and 
equations (10-12), is governed by the initial conditions imposed to the transition by the 
sheath. The ultimate condition that determines the anode fall characteristics, however, is the 
initial field strength imposed upon the region: if the field is sufficiently strong, the ion density 
is pushed past its well toward quasi-neutrality - these matching conditions, for the transition 
equations, determine the rate of density increase and consequently the initial conditions for 
the ambipolar equations. These criteria will then impose the behavior for the ambipolar 
region, whether the manifolds of the plasma boundary stagnation restrict the system 


characteristics, or whether the realm of instability is entered (see Chapter V and Figure (4)). 
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VH. SUMMARY, CONCLUSIONS AND RECOMMENDATIONS 


A. SUMMARY 

Voltage losses at electrode boundaries, surface erosion and sheath effects play an 
important role in plasma devices and must be controlled in designs of practical interest. 
Particularly, thermal arcjets and MPD accelerators deposit between 1596 and 809^ of the 
input power to the anode. This presents not only a severe performance penalty, but also 
complicates the thermal design problems since the heat thus generated must be radiated away 
from the thruster surfaces. [Ref. 7, 9] A natural question arises: what is the behavior found 
in the anode region? Perhaps a restatement of the qualitative picture of the anode, as offered 


by Ingold [Ref. 43: pp.59-63], serves as an answer: 


1. The anode fall voltage is on the order of the ionization potential of the gas since 
electrons must be accelerated to an energy sufficiently high so that an ion current for 
the positive column is provided. 


2. The anode fall region is determined primarily by space-charge and secondly by the 
ionization requirements. 


Moreover, in a comprehensive review of anode-spot phenomena, Miller [Ref. 44] 
suggests that a sudden increase in the voltage noise frequently indicates that an anode foot- 
print could be forming or that a transition into the anode-spot mode could be occurring. 
Further, diis abrupt change in voltage has been associated with a sudden change of ıon density 
in the region near the anode, a density decrease that would leave the local electron space 
charge uncompensated and thus produce a low-conductivity region. The sudden increase jn 
voltage would then reflect this local shortage of ions - an ion depletion region or ion 
starvation region. Even-though the most suitable model for the anode region is still actively 
debated, the general idea of the appearance of an ion deficiency region near the anode as 
triggering the transition into the anode-spot mode has been fairly supported. 

So to expand these qualitative portraits of the anode fall, this work investigated the 


nature ofthe voltage drops in the vicinity of a non-emitting, positive electrode. The selected 
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approach involves non-linear analysis techniques of the continuum governing equations for 
steady-state, isothermal conditions in one dimension, where both ionization and two-body 
recombination processes are considered. The following conclusions and observations are 
offered: 

B. CONCLUSIONS 

Ions, created through electron bombardment of neutral atoms throughout the anode 
fall region, are repelled electrostatically from the positive electrode toward a location of 
quasi-neutrality at the sheath boundary, as well as toward an unstable stagnation region near 
the plasma interface. These equilibria form as a result of the chemistry present. That is, 
continuity (equations (9)) seems to require flow equilibrium when quasi-neutral conditions 
are initially met within the sheath and when the plasma boundary is reached. It is 
recombination that establishes the latter while the presence of ionization stipulates the former. 
Fundamentally, the sheath, as formulated, is inherently unstable - ions are driven out of the 
system, toward the negative electrode. The stagnation regions were obtained analytically, 
through non-linear analysis techniques; to test these results, nitrogen data were used in a 
numerical algorithm from which the following observations are made: 

At the electrode wall, in its immediate vicinity, at least two major processes repel the 
positive ions away from the wall in the isothermal, steady-state case: the initial electric field 
strength and the initial density gradient. If the field strength is low enough, the gradient 
dominates the system and the ion density behaves as expected, increasing monotonically for 
a short distance. Ifthe field strength is high enough, the gradient influence becomes secondary 
in nature, the field then driving the species densities. However, the low initial field strength 
and ion density gradient are insufficient to push the species toward quasi-neutrality, when 
nn, On the other hand, when the initial field becomes sufficiently strong as compared to 


the species density gradient, the ions are pushed toward what appears as a density trough. 


Furthermore, for certain initial values of An=|n, -n | vs E(0) , the ions are pushed past 


e! 


their well, eventually approaching the electron density in magnitude but not merging: for 


sufficiently strong initial field, quasi-neutrality forms in a self-consistent way. Further, if 
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conditions permit the 10ns to traverse their density well, their constriction, then at the moment 


ofquasi-neutrality, when — 7-7, ,there appears a conservative phenomenon for the ions, 


where energy growth is matched by energy loss - a chemical balance. 

At the onset of quasi-neutrality, however, the ensuing species densities do not match 
the plasma neutrality requirements (Saha). Thus, a transition region must drive the system 
toward the plasma boundary: for steady-state conditions, the primary mechanism is the 
density gradient. When the species density becomes large enough, mutual Coulomb fields 
affect free diffusion and ambipolar diffusion results. The ambipolar field then continues to 
drive the ions away from the positive plate; the initial conditions to the transition region, the 
consequent matching to the ambipolar zone, determine the characteristics of this latter region. 
Specifically, the behavior of the anode fall region, as modeled by equations (7-9) and 
equations (10-12), is governed by the initial conditions imposed to the transition by the 
sheath. The ultimate condition that determines the anode fall characteristics, however, is the 
initial field strength imposed upon the region: if the field is sufficiently strong, the ion density 
is pushed past its well toward quasi-neutrality - these matching conditions for the transition 
equations determine the rate of density increase and consequently the initial conditions for the 
ambipolar equations. These criteria then determine whether the manifolds of the plasma 
boundary stagnation restrict the system, or whether the realm of instability is entered. As 
such, any numeric solutions will most likely be challenging, a result observed when, using 
nitrogen data, the various anode fall regions are numerically matched. 

Both the analysis and the numerical results suggest that the anode fall model used in 
this work does not represent sufficient dissipation, preventing the system from adjusting itself, 
to adhere, to the “stable manifold" of the stagnation region near the plasma boundary. The 
model is most likely incomplete in its physical representation: either temperature gradients, 
diffusion written in terms of the electric field, full time dependencies or quasi-steady 
formulation, magnetic field effects and/or three-body recombination effects may provide the 


necessary dissipation. 
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C. RECOMMENDATIONS 
To improve upon the results obtained in this work, further considerations should be 


given to: 


1. Species temperature variation throughout the fall region, specifically in the sheath. 


2. Diffusion perhaps written as a function of the electric field instead of as assumed 
constant for each species. 


3. Full time-dependent and/or quasi-steady formulations for the anode fall. 


4. The effect of three-body recombination on the stagnation regions found, as well 
as the influences of magnetic fields on the system behavior. 


This work presents a broad overview of many subjects, many of which really need 


further investigation, among them: 


1. Analytically matching the second order sheath equations to the transition set and 
to the ambipolar equation; essentially, the analytic analog of the numeric matching 
ideas presented in this work. The results of such mathematical treatment might reveal 
that, in fact, more physics is needed in equations (7-9). 


2. Further investigation of the sensitivity to initial conditions for the ion density well 
as well as quasi-neutral zone formation; a wider range of initial conditions should be 
tested in the hopes of finding some better relation between An versus E which would 
be of aid in any numeric endeavor. Such analysis should also focus on a variety of 
gases, not be limited to nitrogen; possible candidates include xenon and argon. 


3. Generation of a current versus voltage curve, I-V curve, by integrating the electric 


field data and using the current calculations given in this work. Does the result match 
the literature? l 
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APPENDIX A. NON-DIMENSIONALIZATION 


Mathematical derivations are frequently burdened by complicated functions of the 
constants in the problems. Before the mathematical analysis is carried out, a preliminary 
dimensional analysis may not only reveal the ways in which some of the constants enter into 
the final solution, but also unveil significant dimensionless products of a problem. Then the 
original differential equations may be expressed in terms of dimensionless notation, a method 
particularly useful since model laws are usually derived from the differential equations that 
govern phenomena. [Ref. 59: p. 144] Dimensional analysis 1s a process in as much influenced 
by art as it is by scientific methods: two principles were applied to this work - the 
Buckingham Pi Theorem and Fractional Analysis coupled with the physical requirements 
brought by the governing equations themselves. 

A. BUCKINGHAM PI THEOREM 

The Pi Theorem may be stated as “The number of dimensionless products in a 
complete set is equal to the total number of variables minus the rank of their dimensional 
matrix.” [Ref. 59: p. 31] In using the Pi Theorem, the following conditions should be fulfilled: 
[Ref. 60: p. 20] 


1. The list of dimensional parameters must contain all of the parameters of physical 
significance including all independent parameters and one dependent parameter. 


2. The non-dimensional pi's as finally composed should contain, at least once, each 
of the parameters in the original list. 


3. The list of dimensions used to compose the physical parameters must be 
independent, or else provision must be made to compensate for the redundancy. 


There are infinitely many different complete sets of dimensionless products that can 
be formed from a given set of variables. Insofar as Buckingham's theorem is concerned, any 
such complete set is admissible, with some sets more useful than others. So how may a 


complete set of dimensionless products be most advantageously selected at the outset? 
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1. Arrangement of Variables and Terminology 

As Buckingham has pointed out, the maximum amount of experimental control over 
the dimensionless variables is obtained if the dependent variable does not occur in more than 
one dimensionless product, a product referred to as the “dependent dimensionless variable 
or pi" [Ref. 59]. Using a dimensional matrix as outlined in reference 59, the preceding 
condition will be realized, as nearly as possible, if in the dimensional matrix the first variable 
is the dependent variable, the second variable is that which is easiest to regulate 
experimentally, the third variable is that which is next easiest to regulate experimentally, and 
so on. [Ref. 59: p. 39] Consequently, the variables describing equations (7-9) are arranged 
as follows: E=x,, MX, N¿X3, J¿X4, Y=X5, DR, V=X7, IR, DR, KR, Ton, 0X12> 
€,7X,, and q=x,, where E represents the electric field, n,, characterize the ion and electron 
number densities, respectively, j;, are the ion and electron current densities, respectively, y is 
the geometry coordinate length, D,, characterize the ion and electron diffusivities, 
respectively, and v is the ionization coefficient. Universal constants and those variables held 
constant are entered last in the matrix since these are not deemed "experimentally 
controllable”: kis the Boltzmann constant, a, the two-body recombination coefficient, and 
€, the permittivity of free space with q as the charge. The entries of the matrix are next made: 
the dependent variable, i.e., the first entry in the matrix, is the electric field, since in terms of 
voltage applied, it governs the space charge development. The second and third entries into 
this matrix represent the number density distributions and govern the current development. 
All other variable entries are made arbitrarily, the only requirement being that arrangement 
which gives a linearly independent solution for x,4...x4,. To complete the dimensional matrix, 
each variable and constant 1s written in terms of the fundamental units of mass (M), length 
(L), time (t), charge (Q) and temperature (0); the matrix outlining the coefficients of these 


units for each variable becomes: 
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Table 7: Variables in Terms of Fundamental Units 


There will be 14-5=9 dimensionless pi's (number of variables minus the number of units 
necessary to describe them). In fact the rank of this matrix is five which not only reflects the 
number of fundamental units necessary, but also the five constants of the problem: k, T,, q, 
€, €. Solving for the variables x for i-10..14, in terms of the variables x, for j=1..9 results 
in the homogenous, linear, algebraic equations. By taking this equation set into matrix format 
and solving for each variable x; where i=1..9, the homogeneous, linearly independent non- 


dimensional pi’s can be determined: 


Do 2 ty ar X, *X5, 7X5 73x47 5 PF; 
eA. -3% 3x 584 TXC— X 3 X Si 
"loge soe. X. 7 Xo 
A In - 5x, + 3, I, X 


Xj473x1*6 x, *6x,*49x,-2x, *2x,*6 x, 49 x, «2x, 


v) D: (k) | (T) | (2) | (E) 1 (9) 


pedo fe Jefe de fo Tp [> Lp pp 































Table 8: Buckingham Pi's - Matrix Formulation 
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Non-dimensional Relationship Between Non-dimensional | Buckingham Pi's 
Variables Variables, Buckingham Pi's and 
Undisturbed Plasma Values 


N 


L are characteristic 
length scales 


n, = — A 
9 3.3 
(KT, ) "e, e 





Table 9: Relationship Between Non-dimensional Variables, Buckingham Pi's and 
Undisturbed Plasma Values 


When these pi's are inserted into the appropriate governing equations, each term in the 


equation will have a non-dimensional parameter affecting it. For example, the dimensional 


Gauss = -4(n,-n,) written with Buckingham Pi's gives 
MUSS 


12 








or 


T6 dE _ q le 1 q 12 
T dy o^ el. 


Tq Be 
T¿ M¿E 


O 








dy TL, Te, | 





It is the constants and pi's within the brackets that form the non-dimensional parameters. 
2. Pi Theorem - Results and Conclusions 
Although the Pi Theorem forms a good framework for discussing the nature of units, 
dimensions and related topics, the theorem, when used as the only means toward analysis, 


suffers from several deficiencies: 


l. No direct means for finding the pertinent pi's is available and dimensional analysis 
itself provides little framework for incorporating relevant physical information. 


2. The theorem alone does not provide conditions under which one or more pi's can 
be neglected or for which sets of dimensionless groups may be combined. 


So another method must be used to make the problem non-dimensional. 
B. FRACTIONAL ANALYSIS 

There is another common method of dimensional analysis that is intrinsically the same 
as the foregoing technique but differs from it in that the differential equations are expressed 


in dimensionless forms by applying, among others, a characteristic length L and/or a 


RA A — X EX E : 
characteristic time period T, such that zc T-— ;that is, all vanables are non- 
T 


dimensionalized by introducing OAR characteristic values - fractional analysis. [Ref. 
59] A fractional analysis is any procedure for obtaining some information about the answer 
to a problem in the absence of methods or time for finding a complete solution, a method 
from which an approximate analytical or numerical solution can be obtained [Ref. 60]. 


Specifically: 
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1. What is the physical meaning of each of the governing parameters and variable? 
What are the qualitative effects of an increase or decrease in any given parameter or 
variable? 


2. Can we find the conditions under which the effects of certain parameters can be 
neglected either in a given region or for a particular problem? If so, does this lead to 
governing equations that are more tractable so that they can be solved even though 
the general equations cannot be solved? 


3. Are there any combinations of two or more non-dimensional parameters or 
variables, or transformations of variables, which lead to fewer independent quantities 
or which simplify the correlations achieved? 


To answer these questions, the confines of traditional dimensional analysis, as a tool 
by itself, are exceeded and other methods must be employed. One such method incorporates 
dimensional analysis with governing equations as advocated by Sedov [Ref. 61]. Sedov has 
clearly shown that important information can be obtained by simultaneous use of the 
governing equations together with dimensional analysis. Consequently a systematic 
methodology is needed to obtain information directly from the governing equations without 


actually solving them, a methodology involving three primary ideas: [Ref. 60: p. 66-70] 


1. Normalization based on governing equations; making equations and conditions 
non-dimensional in terms of non-dimensional variables of standard magnitude. 


2. Absorption of parameters. 


3. Combination of variables. 


To carry out such procedures, a clear understanding of the physical information 
inherent in the equations is essential. Of equal importance is using a standard procedures for 
transforming the variables to non-dimensional form and standard magnitude; it is not 
sufficient to make governing equations non-dimensional in any arbitrary way. So of the three 
processes stated above, normalization is the most important. [Ref. 60] 


Here, normalization is defined as making the governing equations and conditions non- 
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dimensional in terms of non-dimensional variables of standard magnitude. To carry out a 
normalization, two steps are usually required: 1) making all the variables non-dimensional in 
terms of the appropriate scales of the problem; 2) dividing through the equation by the 
coefficient of one term to make the equation dimensionless (unit-free) term by term. Further, 
the method for choosing the scaling is critical as this choice will permutate throughout the 


ensuing analysis. To this end, the following procedure is used (after Kline): 


l. Define all dependent non-dimensional variables so that they are approximately 
unity over a finite distance and nowhere exceed approximately unity in the domain of 
concern; in this work, the undisturbed plasma values are used, denoted by a “o” 
subscript. 


2. Define all independent non-dimensional variables so that their increment is 
approximately unity over the same domain of concern (0 to 1, 1 to 2, etc, in the new 
variables). 


As a result, the dimensionless groups formed are usually composed from the boundary 
conditions; from the charactenzing sizes or scales of the body; and from the physical 
constraints of the original equation such as system properties, physical constants or both. 
In general, it is not known in advance, which choice of pi's gives the most useful set of 
parameters but it is obvious that the fewer pi's required to specify the behavior, the more 
useful the result will be. So, how to reduce the number of pi’s and consequently parameters? 
[Ref. 60: pp. 75-94] The choice of length or time-scales coupled with physical insight are the 
first step in this effort, an effort leading to a “modified fractional analysis.” 

C. MODIFIED FRACTIONAL ANALYSIS: ANODE FALL 

]. Non-dimensional Parameters 

The derivative term found in the governing equations for the anode fall, equations 
(7-9), is important since without which the conditions at the plasma boundary could not be 
satisfied. So to make all the variables non-dimensional in terms of the appropriate scale, the 
equations are divided through by the coefficients of the derivative terms to form the necessary 
system parameters where both Buckingham Pi's and undisturbed plasma values were used to 


normalize the variables. Following is a list of the results obtained: 
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Governing Equations Non-dimensional Non-dimensional 
[Equations (7-9)] Parameters via Parameters via Undisturbed 
Plasma Values 


where Die = D; e 


ITA n, L,| 


j= = - | 





Non-dimensional Parameters 


2: Estimate of the Length Scales 

In the above, the variables where normalized using undisturbed plasma values since 
the Buckingham's Theorem does not reflect the physics at hand: for example, the process of 
ionization affects the number density distribution in the sheath resulting in electric field 
variations and consequently influencing the current distribution. It would be difficult for one 
length scale to capture each of these behaviors. Consequently physical reasoning must be the 
foundation from which to estimate the length scales of the problem. 

The neutral plasma is shielded from the space-charge region through coulombic 
effects, effects which should drive the length scales affecting the associated momentum 
exchange processes. The first term on the right-hand side of equations (7a) and (7b) model 


these effects and are used to estimate the length scale affecting the number density 


distribution:  L,,L,  . Specifically, 
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E =|a|nE- || D. Equation (7a) 





d(n,) --[«] zz «[a je Equation (7b) 


where [a] and [c] are the coefficients for the coulombic terms summarized in Table 2 (main 


body of this work). Through the process of non-dimensionalization, the parameters are O(1) 


kT 
and the length scales become: Z, As e = . Similarly, the length scale for the 
' q q 


o 





o 


electric field comes strictly from Gauss, equation (8), where either term can be used to 
estimate the electric field length scale since the undisturbed plasma value for n, 1s the same 


for both 1ons and electrons: 


= = E ln, -| fJ, Equation (8) 


^ 





: : ae. e E 
Using the non-dimensionalization result that all parameters are O(1), — L,-——* ^ where 
qn 


M 


the parameters [e] and [f] are again given in Table 2. 
Within the sheath, the ionization process outweighs the recombination loss, so in 
equations (9a) and (9b) the ionization term, the first term on the right-hand side of the 


equations, is used to estimate the current density length scales: 





mee! -|e |v, +% 7,7, Equation (9a) 


202.[ s; jg, ^ Equation (9b) 
That is, the non-dimensional parameters [g] and [i] (Table 2) estimate the length scales for 


10 


the appropriate current density: L ——C-—, L -——-—  . Essentially, the current scale 


reflects the influence of ionization. 

3. Dimensionless Products - Groupings 

Occasionally transforming dimensionless products to achieve greater control of the 
variables is desired or may form a useful grouping: the various results formed through non- 
dimensionalizing the Navier-Stokes equations are examples. The following combination of 


parameters have been found from equations (10-12): 
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Parameter Groupings 


from Equations (10-12) 






Groupings in Terms of Problem Constants 





Q 
Q 
nn; 





Table 10: Non-dimensional Parameter Groupings 
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4. Summary 

Fractional analysis coupled with the use of the governing equations indicates that 
potentially five length scales are present, one for each of the governing equations; however, 
by using undisturbed plasma values for the normalization, three length scales result: a scale 
reflecting the conservation of momentum, one for the variation of the electric field, and the 
third for charge conservation. Of these three lengths, the electric field scale is the smallest, 
indicating the extent of the sheath as the space-charge region defines this region. The number 
density scale, on the other hand, acts over the entire range of the anode fall and so should be 
larger in magnitude than that for the electric field scale. The largest of the length scales is the 
current scale since the current follows the variation in species density according to j,7qn,v, 


where s-i,e: 


Sheath: Transition Zone: Ambipolar Region: 


Space-Charge Dominant Number Density Gradient Equilibrium Near Plasma 


Dominant 





Length Scales via Fractional Analysis 


80 


APPENDIX B. IONIZATION AND TWO-BODY RECOMBINATION DATA 
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= Figure 27: Two-Body Recombination 
Figure 26: Ionization Coefficient v/N as Coefficient as Function of E/N for N, 
Function of E/N for N, [after Ref. 41] [Ref. 40] 


Note: 1 Td=10"* [Vm?] 
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APPENDIX C. NON-LINEAR DYNAMICS - AN OVERVIEW 


When exact closed-form solutions of a physical system as modeled by appropriate 
differential equations are difficult to obtain, or when the exact solution is too complicated to 
be useful, then the first step toward an approximate solution is local analysis. The results of 
such analysis are valid only in a sufficiently small neighborhood of a point, so ultimately a 
uniform approximation to the behavior of the solution over an entire interval may be found 
by piecing together regions where the local behavior is known through global analysis 
techniques involving perturbative and asymptotic methods. The following is a brief overview 
and adaptation of references 55 and 54 as well as class notes, reference 62. 

A. LOCAL ANALYSIS: THE BASICS 

If the general system given by 


X =f, (%} x.) 


X, =f (x, uc) 


is visualized as trajectories flowing through an n-dimensional phase or velocity space with 
coordinates (x,,..,x,), differential equations can be interpreted geometrically as vector fields, 
where now the velocity of the flow, the rate of change of the system behavior, can be plotted 
and the corresponding equilibria (stagnation or fixed points) of the system determined where 
the flow velocity is zero: the fixed points are calculated when the first variation tends to 
zero - calculus of variations. Linearization techniques, local Taylor series expansion about 
the stagnation, can now be used to find the stability about these points. Such analysis 
contains qualitative information about the dynamical behavior of the system leading to an 


overall explanation as to “why things happen the way they do.” 
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In general, the behavior of a one-dimensional dynamical system x=f(x) is 


considered as a fluid flowing along the real line with local velocity f(x). This imaginary fluid 


is referred to as the phase fluid and the real line is the phase space. In higher dimensions, a 


similar situation occurs. The solution to — X,zf(x,..,x,) V i-1,2,-,»  givenan 


arbitrary initial condition x, is then represented by the trajectory of the flow in n-dimensional 
space as determined by the function under investigation. The collection of these trajectones 
are referred to as a phase portrait whose appearance is controlled by the fixed points 
(stagnation, equilibria or zeros) x* as defined by /(x*)=0. This fixed point is defined as stable 
if all sufficiently small disturbances away from it dampen out as t>~. Conversely, unstable 
equilibria are those for which disturbances grow in time. 
B. LINEAR STABILITY ANALYSIS 

An n-th order differential equation can be expressed as a system of 7 first-order 
equations by appropriate changes in the variables. To determine the rate of decay to a stable 
fixed point, linearization techniques about the equilibria are applied. For example, let x* be 
a fixed point and let n(t)=x(t)-x* be a small perturbation away from x*. To see whether the 


perturbation grows or decays, a differential equation for n is formed: 


._d = HN. $2 i ; 
Tj = ae )=x since x* remains constant by definition of stagnation points. 
í 


Consequently, 


12x" )exefG) mfi n) where fixen) nfa) O) 
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Noting that /(x*)-0, | i12 nf/(x')-*O(m?)  . Moreover, if f (x*)«0, then 


wn n is the linearization about x* where the perturbation n(t) grows 


exponentially if /’(x*)>0 and decays exponentially provided / ’(x*)<O. In essence, the slope 
or the rate of change of f '(x*) at the fixed point determines the stability of the stagnation. 
If the dimension of the dynamical system is more than one dimensional, the procedure 
for stability analysis is similar but now a multi-variable Taylor series expansion is used for the 
linearization technique. Consider a two-dimensional non-linear system where the general 


form of a vector field on the phase plane is given by 


ži =J px) 
DUM) 


with f, and f; as given functions. In vector notation, x=(x, ,x,) and / (x) 7 [f, (x), (x)]. Here, 


x represents a point in the phase plane and X is the velocity vector at that point. By 


flowing along the vector field, a phase point traces out a solution x(t), corresponding to a 
trajectory winding through the phase plane. Essentially the entire phase plane ıs filled with 
trajectories since each point can play the role of an initial condition. For nonlinear systems, 
however, it remains a challenging task to find these trajectories, the solutions, analytically. 
Even when explicit solutions are available, their form is often intractable and provide limited 
insight: the importance of a qualitative analysis to understand the characteristics, the 
dynamical behavior, of a system becomes apparent. 


Consider the system 


x= f(x,y) 
Y=g(x,y) 
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and suppose that (x*,y*) represents an equilibrium position, 1.e., f (x*,y*)-0 and g (x*,y*)-0. 
To begin linearization, let u=x-x* and let w=y-y* be components of small perturbations away 
from the stagnation. To see whether the disturbance grows or decays, a differential equation 
for u and w must be found similarly to the one-dimensional method above (noting that x* and 


y* remain constant by definition of equilibrium): 


u=x=f(x*+u,y *+w) 
fe y") up) + (9)+0(u2,w2uw) 
Ox Qy 


EG zn +O(u?,w?,uw) 
ox oy 


and 
W-yog(x «uy ew) 
-g(x' y") +4 (g)+w (g)+0(u2,w2 uw) 
ox oy 


=uL(g)+w (g)+0(u?,w2 uw) 
ox oy 


where partial derivatives 1n the preceding are all evaluated at (x*,y*). That is, the effects of 
the perturbations (u,w) evolve according to the linearized system: 

O Q 
" mU, ay | 
Ww 





u 

| | * O(u?,w?, uw) 
2 (g) 2 (g) m 
ex Qy (x.y ) 

The matrix of partial derivatives evaluated at the stagnation (x*,y*) is termed the Jacobian 
matrix, the multi-variable analogue of the derivative f *(x*) developed for the one-dimensional 
system. Furthermore, if the quadratic terms found in O(1£Z, w^, uw) are sufficiently small, then 
the behavior of the linearized system can be analyzed through study of the Jacobian matrix 


alone. 
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C: CLASSIFICATIONS OF FIXED POINTS 

To classify all possible phase portraits that can occur given an arbitrary Jacobian 
matrix, it should be noted that the coordinate axes play a critical geometric role; they 
determine the direction of the flow as t~+, containing straight-line trajectories. A trajectory 
starting on one of these axes stays on that axis forever and exhibits simple exponential growth 


or decay along it. For the general case, however, the analog of these straight-line trajectories 


has to be found, so trajectories of the form x(1)-e^'V  , where V#0 is some fixed vector, 


are sought with A as a growth rate, also to be determined. If such solutions exist, they 
correspond to exponential motion along the line spanned by the vector V. The conditions on 
A and V are determined by the eigenvalues and corresponding eigenvectors of the Jacobian 
matrix, respectively. That is, the characteristics of the linearized system are determined by 
the nature of the Jacobian eigenvalues. 

de Robustness of Linearization 

The stability of a fixed point can be classified in terms of robust or marginal cases. 
Robust cases are those for which the stability does not change under the influence of small 
non-linear terms such as repellers (sources), attractors (sinks) or saddles. Specifically, 
repellers are those fixed points for which the Jacobian displays positive real parts for all 
eigenvalues, attractors are those characteristics which display negative real parts for all 
eigenvalues. Saddles occur when one eigenvalue is real and positive while the other 1s real 
and negative. Marginal stability results when both eigenvalues are imaginary (centers) or 
when at least one of the eigenvalues is zero. For marginal stability, then, at least one of 
eigenvalues satisfies Re(A)=0 whereas for robust stability, Re(A)#0; that is, for marginal 
stability, the trace of the Jacobian matrix is zero or the flow divergence is zero, implying that 
the flow volume is preserved. 

2. Real, Distinct Eigenvalues 

In general, trajectories approach the origin of the phase portrait tangent to the slow 
eigendirection, the direction spanned by the eigenvector with the smaller |A|. If |A,]<]A,| , then 


the trajectories approach A, as t-«: 
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Real, Distinct Eigenvalues - 
0<|A,|<|A,| and A, A,<0 
“Node” 


If all eigenvalues are positive and distinct, then the behavior locally about the respective 
stagnation is termed “unstable node”, whereas if all eigenvalues are negative and distinct, then 
the dynamics about the fixed point is called “stable node". 
What happens when one of the eigenvalues is positive while the other is negative? 

If A70, the corresponding eigensolution grows exponentially; if A<0, the eigensolution decays 
exponentially. The stable manifold corresponds to the line spanned by the eigenvector formed 
by A<0, whereas the unstable manifold corresponds to the line spanned by the eigenvector 
formed by A>0. The combination of such behavior forms the “saddle” where trajectories are 


attracted toward the stable manifold and repelled as the unstable manifold is approached: 





Real, Distinct Eigenvalues - 
A,<0, 4,70 “Saddle” 
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3. Complex Eigenvalues 


The eigenvalues of the Jacobian determine the local dynamical behavior of the system: 


the eigensolution x(f)=c,e d V FOE ^ V, is affected by the eigenvalues according 


to the behavior of e" . If A, „=a+iw, then the system is influenced by e“e® where a is the real 
component of the complex eigenvalue (damping) and w is the imaginary component, #0. 
By Euler’s formula, e=cos(wt)+i[sin(wt)] and so x(t) becomes a cosinusoidal function with 
exponential growth, instability, according to e" if a>0, or with exponential decay, stability, 
if «0. Such behavior is termed "spiral", unstable or stable depending on whether energy is 


supplied to the flow or taken from the flow. 
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Complex Eigenvalues: 
A, =0-ti0t “Spiral” 


If, on the other hand, the eigenvalues are purely imaginary, «—0, then all solutions are 


periodic with period T=2r/w where the oscillations have constant amplitude. 





Imaginary Eigenvalues: 
A, 7 =+iwt “Center” 


Centers, or systems with purely imaginary eigenvalues, occur generally where energy 
Is conserved, corresponding to a state of neutrally stable equilibrium, or a position of 
minimum energy, although the linearization represented is not robust. Small non-linear terms 
can change the characteristics from periodicity to exponential growth or decay (see below). 

4. Real, Equal Eigenvalues 

Suppose that A,=/,=A, the eigenvalues of the system are equal. Two possibilities now 
exist: either there are two independent manifolds corresponding to A, or there is only one. 

If two independent manifolds exist, then they span the plane and so every vector forms 
a manifold with this same eigenvalue A. If 1.0, all trajectories are straight lines through the 
corresponding equilibrium position which is now termed a “star”, stable or unstable if A<O 


or if A70, respectively: 


Real, Equal Eigenvalues: 
A,=A,=A>0 “Star” 


If the eigenspace corresponding to the eigenvalue A is one dimensional, i.e., there is 
only one manifold corresponding to A, the corresponding fixed point is a “degenerate node”: 
as t2 where all trajectories become parallel to the one available eigendirection. Essentially, 
this behavior results when two distinct manifolds are scissored together with some of the 
trajectories becoming trapped in the collapsing region between the two eigendirections, while 
the surviving trajectories are pulled around to form the degenerate node with stability 
determined according to A<O or A>0. The degenerate node is a borderline case between a 
spiral and a node; the trajectories are trying to wind around in a spiral (complex eigenvalues), 


but insufficient energy is present to complete the winding: 
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Real, Equal Eigenvalues with One 
Manifold: A,,=A 


5. Miscellaneous 

a. Small Non-linear Terms Change Linearization Stability-Sometimes 

In the presence of non-linear terms, the character of marginally stable fixed 
points can change: centers, stars and degenerate nodes all satisfy the criterion Re(A)=0. For 
example, if the linearization predicts a center, the influence of non-linear effects causes the 
phase portrait to show a stable or unstable spiral, a stable or unstable node; the presence of 
non-linear terms violate Re(A)=0 by either adding or removing energy from the system, 
driving it away from the conserved state (center). Similarly, stars and degenerate nodes can 
also change their character, but unlike centers, non-linear effects change only the character 
of the star or degenerate node, not their stability: stable star to stable spiral or stable node, 
for example. 

The phase portrait is structurally stable if its topology is not changed by arbitrary small 
perturbations to the vector field: the phase portrait of a saddle point is structurally stable, but 
that of a center is not since an arbitrary small amount of damping converts the center to a 
spiral. Essentially, if the phase portrait changes its topological structure as a parameter is 
varied, a bifurcation is in progress. 

b. Example: A Glider 
Consider a glider flying at speed v at an angle 6 to the horizontal. Its motion 


is governed approximately by the dimensionless equations 


9] 


v . .sin(8) -Dv? 
dt 


v 98. cos(0) «v? 
dt 


where the trigonometric terms represent the effects of gravity and the v^ terms represent the 
effects of drag and lift. Here, it is desired to find the effects on the glider performance as the 
drag parameter, D, is varied. 

To begin the analysis, first a local analysis is performed about any stationary points 
(fixed points) that occur when the derivative terms are set to zero. These fixed points are 
then used to evaluate the linearization (Jacobian) so that the system characteristics about the 
stagnation can be determined in terms of eigenvalues and eigenvectors. The results are 


summarized: 


Jacobian Eigenvalues 


n=0,2,4,6,... 
Az 2D+/4D?-8 


=> 
172 
j 2 


j .2D*/4D?-8 


12 
2 





From the preceding, it is seen that if the drag D>2 in non-dimensional terms, the eigenvalues 
of both fixed points will always be real and distinct. The associated dynamical behavior, 
locally, can therefore be characterized as a stable node for the fixed point (n7,1) and as an 
unstable node for (nz,-1). If D20, then the trace of the Jacobian (sum of respective 
eigenvalues) will always be less than zero for the fixed point (n7,1), while the instability for 


(n7,-1) is shown by the exponential growth of the oscillation component of the associated 
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system eigenvalues (trace>0). 

If 2>D20, the eigenvalues are complex and hence the local dynamics of the system 
is represented by stable spirals for (nz,1) and an unstable spiral for (nm,-1). Similarly, the 
stability/instability is evident through the exponential decay/growth of the oscillations 
characterized by the real component of the respective eigenvalues. It should be noted that 
when D=2, the eigenvalues are real and equal, however. Consequently, the system behavior 
should be that of a star or a degenerate node: a spiral (unstable or stable, depending on the 
eigenvalues), 1s found instead. Here is evidence that nonlinear terms in the system can change 
the local behavior from that which is predicted by the linearization analysis. 

If D=0, the case of no drag, then the eigenvalues for the fixed point (nr,1) with n=0 
will be imaginary, reflecting pure oscillation of the system and consequently periodic solutions 
with period T-27/0. The resulting equilibria are characterized as centers and describe the 
lowest energy or power state of the system, the glider moving in a potential field, gravity, 
without any other influences. That is, the glider moves along lines of constant energy, along 
a curved flight path with a radius of curvature equal to the absolute altitude, experiencing the 
standard ^ fictitious" forces due to relative motion (centripetal acceleration, for example). 

In general, for D>0, the case of increasing drag, the centers at (n7,1) transition to 
stable spirals/nodes, whereas the centers at (nrt,-1) transition to unstable spirals/nodes. From 
an eigenvalue consideration, both dampening and oscillations are now evident; the eigenvalues 
are complex with exponential decay for the stable spiral/node and exponential growth for the 
unstable spiral/node. In essence, drag tends to offset the lift generated, as expected. This 
example shows that if the phase portrait changes its topological structure as a parameter is 
varied, a bifurcation 1s in progress. 

D. © BIFURCATIONS 

As indicated, the stability of a fixed point can change; such qualitative changes in 
dynamics are referred to as bifurcations. Bifurcations provide models of transitions and 
instabilities as some control parameter is varied, the onset of coherent radiation in a laser, for 
example. Of the many types of bifurcations possible, only Saddle-node and Transcritical, are 


outlined here. 
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The saddle and transcritical bifurcations occur when one of the eigenvalues equals 
zero, involving the collision of two or more fixed points. These types of bifurcations have 
analogs in one, two and higher dimensions since the important behavior for these dynamics 
are confined to a one-dimensional subspace along which the bifurcation occurs, when 
Re(A)=0, while in the extra dimensions the flow is either simple attraction or repulsion from 
that subspace. 

fh Saddle-Node Bifurcations 

The saddle-node bifurcation is the basic mechanism by which fixed points are created 
and destroyed. For instance, as a parameter is varied, two fixed points move toward each 


other, collide and mutually annthilate. The prototypical example of a saddle-node bifurcation 


X=r+x? 


E where r is a small parameter determined by the 


is given by the system 


problem at hand. Here, the motion is decoupled with the y-direction assumed arbitrarily as 


exponentially decaying. The fixed points for this system are — (x y )- (er) , 


stagnation points which exist in the x,y plane when r>0, coalesce when r=0 and do not exist 
(in x,y space) when r<0. Upon further analysis, it is seen that the bifurcation is fundamentally 
one-dimensional, with the fixed points sliding toward each other along the unstable manifoid 


ofthe saddle point at (x ’,y)=(-Yr,0) where r>0. 


For two-dimensional system, and for higher-dimensional systems, the flow is limited 
not only by fixed points (stagnation), buts also by closed orbits and the unions of fixed points 


and the trajectories connecting them. The latter are referred to as heteroclinic orbits when 
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they connect distinct points and homoclinic orbits when they connect a point to itself. For 
example, if the physical system is such that the unstable manifold of a saddle node intersects 
with its stable manifold, the resulting orbit ıs called “homoclinic:” the connection of the 
unstable manifold (A>0) and the stable manifold (A<0) of the saddle point is common in 
conservative systems but can also occur in non-conservative systems. When the system 
dynamics are such that additional damping is introduced, the tendency is for the homoclinic 
orbit to break, causing the unstable manifold to have components which approach the fixed 
point at the origin of the phase portrait as t>~. This fixed point is then a sink, with 
eigenvalues -a+i\ where a is the damping. [Ref. 55: p. 45] 

2: Transcritical Bifurcation 

There are certain scientific situations where a fixed point (stagnation) must exist for 
all values of a parameter and can never be destroyed. For example, in the logistic equation 
there is a fixed point at zero population, regardless of the value of the growth rate; the laser 
rate equation (classical) is another example for which a fixed point always exists. However, 
such a fixed point may change its stability as the parameter is varied. The transcritical 


bifurcation is the standard mechanism for such changes in stability, with normal form given 


oo 4 Pe 
by: xX-FXIX 


is where r is control parameter determined by the problem. The 


difference between the saddle-node and transcritical bifurcations is that in the latter, the fixed 
points do not “disappear,” their stability changes instead. Consider, for example, the laser 
rate equation [Ref. 54: p. 54]: 
= -gain - loss- GnW-kn 
í 


=Gn(W_-an)-kn= GW n - aGn° -kn 
=(GW_-k)n-(aG)n* 


where n(t) represents the number of photons, W characterizes the rate with which a single 


excited atom spontaneously generates a photon per second offset by the number of photons 


25 


per second that drop back to the ground state, a; W, 1s the pump strength. G is the gain 
coefficient, and -K shows the rate at which photons are lost in the laser by scattering or 
impurities. For the laser rate equation. shown, two equilibria exist: n,*=0 and 
. GW,-k dine 
n, = — But what does this mean? Shown below are the phase portraits with 

c 


the fixed points indicated: 
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When W,< k/G, GW,-k<0, n,*=0 is stable and for low levels of incident energy (weak 
pumping) the solution n,*=0 must be valid since n(t), the number of photons generated, must 
always be positive. Consequently, there can be no laser action - the laser acts like a lamp. 
When the laser is pumped at a greater intensity such that now (GW,-k)>0 or W,>k/G holds, 
the solution with n>0 becomes possible and laser action is attained. The system undergoes 
a transcritical bifurcation when W,=k/G, the laser threshold for this model. When W,>k/G 
the lamp becomes a laser and the fixed point n,*=0 loses stability, driving the system toward 
n,*. The pumping intensity elevates the inversion population: a critical threshold exists when 
W,-k/G below which laser action cannot occur (insufficient pumping energy and 
consequently insufficient inversion) and the laser acts like a lamp. Above the threshold, laser 


action occurs since sufficient energy is now available to raise the atoms to higher energy 
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levels so that the number of photons generated are enough to cause lasing. At the threshold, 
at the transcritical bifurcation point, the stability of the system origin becomes unstable, 
driving the dynamics toward greater population inversion and consequently toward the second 
fixed point. In fact, there are certain scientific situations where a fixed point (stagnation) must 
exist for all values of a parameter and can never be destroyed although it changes its stability 
character; the laser is such an example, for it is the stagnation at the origin that drives the 


system toward the lasing threshold at the bifurcation point. [Ref. 63] 


o" 





APPENDIX D. CHARACTERISTICS OF SHEATH EQUATIONS 


The Jacobian to equations (7-9), 


= — b 

È | = 
[a] 0 [a]n, 0 EA 

area ae 

0 [cla [c]n, M 0 

Å inner = [e] - 1 0 0 0 

(Ala, -[glv+[aya, ED 0 

us 
A O 
-Uln, Ulv-pl] ———— 0 


r2 


has the following characteristic polynomial: 


characteristic polynomial-aX + BA + yA? + 212 + n2+€=0 


where 


99 


a=D,D,E* Pß=D,D,E’\c-a] 
v-E'D- [an] A P E 
E? D, |»; |, + E'D,D, - ea] - GA = 
E*D,Dac| 
C-E^Djbjc n, + E’ D aana : lade] : 
E'D,D,l-|eac]r, + lare]. + 
AD, ebilv - n,D|fed|v 
n=En,D,ecbilv + 
Es, | [anbi - [aes] + 
En,D,|afdg |v + 
72 E*D,| bjef E 
E^D nx |anaf - n|aedg v + à [acean] + 


E'D, i [reb - 2 


e- i [ebjgali : ledhiblv - n hdfib v - [areas 


To find the solutions to the characteristic polynomial, at least one of the roots has to 


be known apriori. To this end, the Jacobian is evaluated for the fixed point found when 


n, =n*=n"=0 (see Table 5 in the main body of this work) but with E retained as 


e 


a non-zero variable (see Chapter V). The resulting Jacobian eigenvalues for the fixed point 


for this condition become: 
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To find the general eigenvalues of the Jacobian then, the characteristics 


A,-[a]E and A420 are assumed as solutions. By factoring the characteristic 


polynomial with these results and neglecting any remainder terms gives the approximate and 


complicated solutions to the characteristic polynomial: 








A, =[a E 2,=0 
Jy3 +x2 > ee 
SE Y P Y poz 
x x 
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3\ a 
where 
x-|a E, v- +x) x + =, Y =YıK + £ 

a a 


If numeric data nitrogen is applied to the stagnation region near the plasma boundary (Tables 
4, 5) and the results substituted into these roots, an order of magnitude comparison can be 
performed with the actual results obtained (the exact evaluation of Jacobian at the fixed point 


in question): 


E EE 


1.376704 » 0.630455 | -1.37437 | -0.63279 
1.41467 0.5 -1.41262 | -0.54205 
to Jacobian at Plasma l 
Boundary Stagnation 


Table 11: Sheath Equations - Jacobian Characteristics Evaluated at Plasma Boundary 
Equilibrium 











Jacobian Evaluated at 
Plasma Boundary 
Stagnation 















Approximate Solutions 





It seems that A, dominates the exponential growth seen during any numeric simulation - the 


process of ion repulsion from the positive plate forces the instability of the sheath. 
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